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Abstract 

A  fine  attituda  dataraination  aodal  is  davalopad  for  a  spinning 
gaosynchronous  satellite,  based  upon  stellar  observations  froa 
a  V-slit  star  scanner  and  ground  supplied  satellite 
epheaerides.  A  true  attituda  aodal  is  datarainad  through 
nuaarical  integration  of  Euler’s  aoaent  aquations  for  a  torque- 
free,  rigid,  axisyaaetric  body,  and  kineaatic  relations  which 
consist  of  pitch,  roll,  yaw,  orientation  angles.  Next,  the 
closed  fora  solutions  to  the  Euler  aquations  are  coupled  with 
first-order  approxiaate  solutions  of  the  kineaatic  relations  to 
develop  a  second  order  kineaatics  aodal.  Observation 
relations,  relating  stellar  slit-plane  crossing  tiaes,  (a 
priori  star  identification  assuaed ) ,  to  attitude  states,  are 
then  developed.  Finally,  a  nonlinear  least-squares  estiaation 
algoritha  is  used  to  identify  the  full  satellite  attitude 
state.  Siaulation  is  tested  for  single  scan  and  aultiple-scan 
capability  using  exact  and  noise  ridden  data.  ?r  i.o:  r\~  ; 


STAR  SENSOR  FINE  ATTITUDE  POINTING  MODEL 
FOR  A  SPINNING  GEOSYNCHRONOUS  SATELLITE 


Chapter  Ona  Introduction 


Topic  Description 

Today* a  aatallitaa  utiliza  a  variaty  of  aenaora  and 
inatruaanta  for  attituda  determination.  For  a  typical 
spinning  aatallita,  data  froa  a  tachoaatar,  aun  sensor, 
aarth  sanaor,  and  atar  aanaor  would  ba  couplad  with  an 
aphaaaria  pradiction  aodal  and  fed  through  a  ground-based 
computer’s  attitude  prediction  algorithm  to  yield  the 
satellite’s  state  ( pitch, ro 1 1 , yaw,  and  rates).  In  the 
quest  for  autonomous  satellite  navigation,  simple  and 
reliable  models  that  accurately  compute  satellite  attitude 
need  to  be  developed.  This  report  analyzes  an  attitude 
prediction  aodal  for  a  torque-free,  spinning,  rigid, 
axisyaaetric,  geosynchronous  satellite  based  solely  on 
satellite  epheaerides,  a  star  catalogue,  and  data  obtained 
froa  a  scanning  star  sensor.  As  onboard  satellite 
computing  capacity  continues  to  increase,  it  is  quite 


feasible  that  satellite-borne  attitude  determination 


schemes  will  soon  replace  ground  based  systeas.  For 
aaxiaua  reliability,  these  algorithas  will  require  as  auch 
independence  as  possible.  Optiaal  pointing  accuracy  will 
be  achieved  through  a  weighted  reliance  on  all  of  the 
sensors,  but  in  the  event  of  satellite  software  or  hardware 
anoaolies,  standalone  capabilities  would  prove  invaluable. 
An  algorithm  siailar  to  the  one  presented  in  this  paper 
aight  one  day  sake  up  an  independent  portion  of  a 
satellite’s  attitude  determination  package. 

Background 

Satellite  star  sensors  were  conceived  to  assist  in 
determining  a  satellite’s  attitude.  Many  of  the  early 
ideas  for  solution  to  the  general  problem  of  spacecraft 
attitude  determination  were  probably  first  consolidated  in 
open  literature  in  the  "Proceedings  of  the  Symposium  on 
Spacecraft  Attitude  Determination"  (It  1-438).  In  that 
report,  attitude  determination  based  on  stellar 
observations  is  a  pervasive  topic.  The  current  need  for 
efficient,  independent  computer  programs  for  possible 
space-borne  application  gives  the  subject  a  renewed 
importance. 

Mode  1  Description 

Any  attitude  determination  scheme  based  on  stellar 
observations  consists  of  an  observer,  a  dynamics  package, 


and  an  estimator.  The  estimator  fits  the  data  recorded  by 


the  observer  with  the  motions  described  by  the  dynamics 
package  and  produces  a  best  estinate  of  the  attitude.  In 
computer  simulations,  the  functions  of  the  observer  Bust  be 
augmented  with  a  truth  model,  which  tracks  the  attitude  of 
the  true  state  of  the  simulated  satellite. 

The  observer  can  be  a  camera,  a  star  tracker,  or  a  star 
scanner.  Star  trackers  operate  by  mechanically  tracking  a 
particular  star.  Data  is  recorded  as  the  orientation 
between  the  tracker  boresight  axis  and  a  satellite  fixed 
reference  frame.  Star  scanners,  on  the  other  hand,  are 
generally  fixed  to  the  satellite,  and  they  scan  the  heavens 
as  the  satellite  moves.  Data  is  recorded  as  the  time  when 
a  star  passes  through  the  f ield-of-view  of  the  scanner. 

The  f ield-of-view  usually  contains  one  or  more  planar 
slits,  called  slit  planes  (2;582).  For  spinning 
satellites,  the  practical  choice  of  star  sensors  is  the 
scanner.  Herein,  a  star  scanner  is  modelled  with  the 
following  specif icationst  +8.0  Magnitude  Sensitivity;  3° 
Field-of-view;  .205’  Accuracy  Ob’);  38°/sec  Scan  Rate. 
Comparing  this  with  a  first  generation  star  sensor,  the 
0S0-8  Star  Scanner,  with  +4.0  Magnitude  Sensitivity;  10° 
Field-of-view;  8’  Accuracy  (3b);  30°/sec  Scan  Rate,  (3s 
224),  it  can  be  seen  that  the  model  specif icat ions  are  more 
stringent,  but  not  unfeasible.  The  scan  rate  mentioned 
above  can  be  either  a  satellite  characteristic,  or,  in  the 


case  where  the  scanner  is  motor-driven,  a  characteristic  of 
the  coabined  satellite  dynamics  and  motor  rate.  A  complete 
description  of  the  star  sensor  model  is  given  later  in  this 
report. 

The  attitude  dynamics  package  describes  the  motion  of 
the  satellite.  A  dynamics  package  may  be  a  piece  of 
hardware  or  a  mathematical  model.  An  example  of  a  coabined 
hardware/software  dynamics  package  would  be  an  aircraft’s 
inertial  navigation  system.  For  this  report,  a 
mathematical  dynamics  package  was  chosen  to  model.  The 
dynamic  equations  consist  of  Euler’s  first  order  moment 
equations  and  kinematic  equations  which  contain  Euler 
orientation  angles.  Many  papers  written  on  this  subject 
use  Euler  parameters  instead  of  Euler  orientation  angles  to 
describe  the  satellite's  kinematics.  While  this  choice 
avoids  singularities  which  are  inherent  in  the  use  of  Euler 
orientation  angles,  it  forces  the  use  of  an  additional 
variable,  and  it  removes  physical  clarity  from  the  dynamic 
description.  This  report  uses  the  familiar  pitch,  roll, 
and  yaw  orientation  angles,  and  while  it  is  acknowledged 
no  solutions  are  possible  for  pitch  angles  of  90°,  later 
assumptions  prohibit  this  situation. 

The  assumptions  of  an  axisymaetric  rigid  body  in 
torque-free  motion  permit  the  solution  of  Euler’s  equations 
in  closed  fora.  It  has  been  shown  that  making  an 


axiayaaatric  aaauaption  in  naar~axiayaaatric  aituationa 
laada  to  vary  aaall  arrora  (2i566).  Tha  kinaaatic 
aquationa  ara  naxt  aolvad  to  firat  ordar  by  iapoaing  aaall 
angla  raatrictiona  on  tha  aatallita’a  pitch,  apin  apaad 
daviation,  praeaaaion  rata,  pracaaaion  angla,  and  orbital 
rata.  For  tha  problaa  of  f ina-attituda  pointing  on  a 
noainally  aarth  pointing  aatallita,  thaaa  aaauaptiona  ara 
not  unraaaonabla.  Tha  firat  ordar  kinaaatic  aolutiona  ara 
than  uaad  to  dariva  aacond  ordar  aolutiona,  which  aaka  up 
tha  aacond  half  of  tha  dynaaica  packaga.  Tha  ganaral 
orbital  dynaaica  problaa  ia  aaauaad  to  ba  aolvad 
indapandant ly  for  thia  aiaulation.  A  full  daacriptlon  of 
tha  attituda  dynaaica  routina  ia  givan  in  tha  naxt  aaction 
of  thia  raport. 

Tha  final  coaponant  of  tha  attituda  dataraination 
ayataa  ia  tha  aatiaator.  Harain,  a  nonlinaar  laaat  aquaraa 
routina  ia  chosan.  It  will  ba  ahown  latar  that  thia 
routina  ia  wall  auitad  for  tha  problaa  at  hand. 

Tha  truth  aodal  waa  aiaulatad  by  nuaarical  intagrationa 
of  tha  axact  dynaaic  aquationa  which  wara  fad  through  tha 
obaarvatlon  ralationa. 


Chapter  Two  Analytical  Developaant 


Ovarviaw 

Tha  coaputar  aodal  for  attitude  dataraination  will 
contain,  a a  atatad  earlier,  four  basic  parte:  tha  dynaaics 
aodal,  tha  observer  aodal,  tha  truth  aodal,  and  the 
astiaation  algorithm. 

Tha  aquations  comprising  tha  dynaaics  aodal  are  Euler’s 
aoaant  aquations,  which  can  be  solvad  exactly  for  tha 
torqua-fraa  rigid  body,  and  a  sat  of  kinematic  relations, 
whose  solutions  are  derived  through  second  order 
approxiaations.  Next,  tha  observation  relations, 
equivalent  to  taking  the  inner  product  of  the  boresight 
axis  of  each  slit  plane  of  the  star  sensor  with  a  star 
vector  represented  in  the  body  frame,  are  derived.  Then, 
the  truth  aodel,  consisting  of  a  bright  star  catalogue  (4< 
H1-H31  ft  SiSec.  VI)  and  a  nuaerical  integration  of  the 
exact  state  equations,  is  fully  described.  Finally,  the 
nonlinear  least  squares  estiaation  algorithm,  which  takes 
Inputs  from  all  other  components  to  produce  a  state 
estimate,  is  explained. 

Satellite  Dynamics 

Coordinate  Systems 


In  order  to  define  the  satellite's  dynaaics,  two 


coordinate  systems  are  employed.  The  first  reference 
fraae,  the  orbital  frame  (0^,02,03),  centered  at  the 
satellite’s  mass  center,  consists  of  a  pair  of  axes  in  the 
orbital  plane  (og  "down”  froa  the  satellite  to  the  earth’s 
center,  and  o^  ’’eastward”  along  the  orbital  path),  and  an 
og  axis  normal  to  the  orbital  plane  (see  Figure  la).  The 
second  reference  fraae  consists  of  a  principal  body  axis 
set  (blfb2«bg).  bg  is  the  symmetry  axis  of  the  satellite 
(nominally  down-pointing),  and  b^  and  b2  complete  a  right- 
handed  set  in  the  spin  plane  of  the  satellite  (see  Figure 


Kinagaticg 


The  kineaatic  equations  of  notion  ara  derived  froa  a 
pitch  roll  (4*i>t  yaw  ((113)  rotation  C2,l,3]  of  the 

orbital  axes.  Thust 


C  " 
*»1 

*11 

*12 

*13 

°1 

t>2 

*21 

*22 

*23 
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02 
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*31 

*32 

*33_ 

°3 

where i 


•ll  »  sin < )sin( 4>2 lain ( *|>3 )  +  cos  (0»2  >008(013)  (2a) 

«12  a  cos(Mii)sin(M>3)  (2b) 

a13  *  sin(0ii  )cos((i>2)sin(<l>3)  -  sin(0i2)coa(0i3>  (2c) 

*21  a  »in(tj>1  )sin(<i»2>cos((i)g)  -  cos(0i2)ain(0>3>  (2d) 

*22  *  008(01!) cos («p3)  (2e) 

*23  =  sin(0»i  )cos(0>2  )°os  (0>3 )  ♦  sin(0i2>»in(0»3)  (2f) 

*31  31  cos  (Oil  )sin(0>2>  (2g) 

*32  *  ”«in<*i)  (2h) 

*33  *  cos  (Oil  )cos(0i2>  (2i) 


(61 423 ) 


Use  of  this  rotation  scheae  instead  of  the  coaaonly 
used  (3,1,3]  Euler  angle  rotation  aoves  the  singularities 


inherent  in  theae  acheaea  out  of  the  area  of  concern.  (7t 
614) 

The  angular  velocity  of  the  body  axes,  with  reapect  to 
the  orbital  axes,  iat 

JSLlb/o  *  [<i>2Coa((j>^  )ain(u>Q)  +  'j>1coa(t|>g )  3  bj  (3a) 

Jll2b/o  *  C«>2co*^l  )coa(«J»g)  -  <i>iain(<J>3 )  ]  b2  (3b) 

JS3b/o  *  C~*>2*1n(«|>i )  +  «i>3  b3]  (3c) 

( 8 1 469 ) 

Now,  for  a  aiaplified  orbital  aodel,  let  the  aatellite 
be  in  an  equatorial,  circular,  prograde  orbit.  Then,  the 
angular  velocity  of  the  orbital  a xea,  with  reapect  to  an 
inertial  geocentric  aet,  iat 

*  -0  Og  (4) 

or,  expreaaed  in  the  body  fraae,  uaing  equation  (l)i 


w1o/I  «  -Ocoa(i|>1  )ain((t>3) 

A 

bl 

(5a) 

w20/*  *  -Ocoa(«|»i  )coa(t|ig) 

A 

b2 

(6b) 

j«»3o/I  ■  0ain(4»i)  bg 

(5c) 

Adding  equationa  (3)  and  (5)  yielda  the  total  angular 


velocity  of  the  aatellite,  expreaaed  in  the  body  fraaes 


W1  “  4ilb/I  *  t~0co«(«l»i)*in{a»3)  +  igcoald^lain^g) 

•  A 

♦  <|>ico*(<pg)  ]  (6a) 

wg  *  W2b/  *  C-nco»(M>i  )co»(4»3>  +  igcoa(4>^  )co*(4>g) 

-  01ain(M>3)]  bg  (8b) 

wg  «  wgb/I  *  COalnC#!*  -  4>2*in<4>i)  +  b3  (6c) 

•  •  • 

Equation  (6)  can  ba  aolvad  for  ,  <i>g,  and  03,  to 
ylald  tha  kinaaatio  aquations  of  aotion  in  f irat-ordar, 
nonlinaar  foras 

*  u>lcoa(<J>g)  -  wg*in(<pg)  (7a) 

$2  *  Mjain^gj/coa^j)  +  wgcoaOfrg^coa^ )  +  0  (7b) 

0>3  *  uii*in(<|>3)tan(«1 )  +  wgco* (013) tan ( Hij  )  ♦  wg  (7c) 

( 6 1 428 ) 


Eular'a  Equation* 

In  tha  principal  body  axia  aat,  Eular'a  aquations  arcs 


(8a) 

(8b) 

(8c) 

whara  A,  B,  *and  C  ara  tha  principal  aoaanta  of  inartia. 


Aw^  +  (C-B  )ugwg  *  M^ 
Bwg  +  (A  — Ow^wg  ■  Mg 
Cwg  +  (B— A)W^Wg  ■  Mg 


and  M^,  Mg.  and  Mg  ara  tha  aoaant  coaponanta. 

Equation  (8),  whan  aolvad  for  wlt  wg,  wg,  yialda  tha 


attltuda  dynamic  aquation*  of  motion  in  firat  ordar  forms 

■  (  ( B— C ) / A )wg(ug  +  Mj_/A  (9a) 
u»g  *  ( ( A—C ) /B )w^wg  ♦  Mg/B  (9b) 
wg  *  ( ( A— B ) /C ) w^wg  +  Mg/C  (9c) 


For  tha  ca am  undar  atudy,  it  will  ba  aaauaad  that  tha 
aatallita  ia  ayaaatric  about  tha  bg  axia.  Thuas  A  *  B. 

Furtharaora,  if  wa  lat  tha  inartia  ratio  ((A-C)/A)  ba 
rapraaantad  by  k,  than  Eular’a  aquation*  raduca  tos 


wj_  =  kwgbig  +  M^/A 
U£  *  “kwj_wg  +  Hg/ A 
«3  *  m3/C 


(10a) 

(10b) 


(10c) 


Torqua-Fraa  Motion 


In  tha  abaanca  of  axtarnal  torquaa  (M1*Mg*Mg=0) , 
Eular’a  Equation*  can  ba  aolvad  in  cloaad  form 


W1  *  w10coaCwg0k(t-t0)]  +  w20ainCwg0k(t-t0)l 
«2  *  -wio*lntwgok(t-t0)  ]  +'  wg0coatwgo*«(t-t0)  ] 
“3  *  w30 


(11a) 

(lib) 


(11c) 


whara  w2q,  and  wg q  arm  all  conatanta  of  tha 


motion. 


The  dynaaic  nodal  to  ba  uaad  in  thia  work  will  involve 


tha  torqua-fraa  solution  Just  developed. 


I  First  Order  Approxi nations  of  the  Kinenatic  Equations 

The  solutions  to  Euler’s  equations  (11)  can  be  used  to 
»  help  integrate  the  kinenatic  equations  (7),  yielding  the 

total  dynanic  description  of  the  satellite’s  attitude  in 
non-differential  forn.  In  order  to  coaplete  the 
integrations,  the  kinenatic  equations  will  first  be 
sinplified.  With  a  noainally  ”down-pointing”  satellite,  d»i 
and  <|>2<  as  well  as  and  d»g,  will  renain  snail.  Also, 

)  with  the  og  and  bg  axes  remaining  nearly  aligned,  and  uig 

will  have  values  on  the  sane  order  as  d»i  and  d>2  (i.e.j 
snail).  A  final  assunption  will  be  that  orbital  rate  (0) 

)  will  be  coaparitively  snail.  Using  first-order  snail  angle 

approxiaations  for  d>i  and  d>2<  equations  (7)  become* 


a  hl^COS  ( d>g  )  “ 

ugsin ( d>3 ) 

(12a) 

dig  *  wisin(d»3)  * 

WgCOS  ( d*3 )  +  0 

(12b) 

d>3  *  w1d>isin(d»3) 

+  wgd>iCOS(d>3)  +  Wg 

(120 

Neglecting  second  order  terns  in  0,  d>i«  d>2*  «i  •  and  uig, 
equation  (12c),  to  first-order,  becoaest 


Integrating t 


$3  *  4>30  +  “30^t“t0^  (14) 

Tha  results  froa  aquation  (14),  whan  inserted  into  the 
right  hand  sides  of  aquations  (12),  will  give  first-order 
approxi nations  for  <ti^,  cjigi  and  dig « 

<Pl  *  w^cosC«l»3o+ui3o(^“^0^  “  w2*^n^30'4'ta,30^“^0^  (16a) 

ii>2  “  “i*in[i|i3Q+“3o(t-to)  ]  +  (‘*2®^n^30'fw30^“^0 ^  +  ^  (16b) 

<P3  *  “30  (16c) 

Substituting  the  solutions  of  Euler’s  equations  (11) 
into  equations  (16)  yieldst 

*>1  *  “ioco®C<k_1>w30(t~fcO>~'*’303 

+  “2o*inC  (k-1  )“3q(  t-t0)-M>3o3  (18a) 

<j>2  *  -“io»inC(k-l)“3o(t-to)-0»3o3 

+  “20co*t(k“l  >“30*  fc-fcO,"'*,303  +  n  (18b) 

<i>Q  *  “30  (18c) 


These  equations  can  easily  be  integrated  to  yield: 


-d<|>io  C  w10cos < O»30  > ““20* 1  n ( *30  > 3 


<19e) 


|  These  second  order  kineaetic  approximations,  coupled 

with  the  solutions  to  Euler’s  equations  (eq.  7),  complete 
the  dynaaic  description  of  the  satellite  attitude  aodel. 

The  solutions  to  the  dynaaic  equations  can  be  coabined 
to  fora  the  attitude  state  vector  {x>,  consisting  of  state 
eleaents  <*»2«  uig,  4 4»3>T.  The  aatrix  whose 

eleaents  are  aade  up  of  the  partial  derivatives  of  the 
state  vector  with  respect  to  the  state  eleaents  is  called 
the  [4]  aatrix.  The  [4]  aatrix  noraally  plays  a  major  rol 


in  an  estiaation  algoritha.  The  [4]  aatrix  is  not  used  in 


Observation  Geoaetry 

Tha  inartial  raferance  fraae  will  be  defined  by  the 
celestial  sphere.  Let  the  origin  of  the  systea  be  at  the 

satellite’s  aass  center.  The  X  axis  will  point  to  the 
first  point  of  Aries  (y),  the  Y  axis  will  point  east  along 

the  celestial  equator,  and  the  Z  axis  will  point  toward  the 
celestial  north  pole.  Since  all  data  will  consist  of 
distant  star  sightings,  parallax  will  be  ignored  (i.e.-the 
reference  fraae  will  be  treated  as  though  the  origin  were 
at  the  earth’s  center).  Also,  since  the  period  of  analysis 
will  reaain  short,  the  aotion  effects  of  the  earth’s  pole 
(precession,  nutation,  wobble)  will  be  ignored.  Star 
inforaation  will  be  stored  as  right  ascension  (a)  and 


declination  (5).  (see  fig.  2) 


recorded  a s  first  slit  crossing  tine  (t^),  and  dwall  tins 
(td).  A  full  dascription  of  tha  optical  geometry  will  ba 
shown  latar  in  this  saction. 

Looking  at  tha  oriantation  of  tha  orbital  fraaa  at  a 
tiaa  (tgc)  whan  tha  satellite  crosses  the  +X  axis,  the 
following  relationship  is  evident  (see  fig.  8): 


tha  orbital  plana  ia  not  inclinad  with  raapact  to  tha 
calaatial  aquator.  In  aatrix  notation,  aquation  (20) 

bacoaaa « 


°l<tg0>  010 
°2<fcgo>f*  0  0-1 
°3<tgo>  -100 


(21) 


Now,  at  aoaa  arbitrary  tiaa  t,  givan  tha  aatallita’a 
right  aacanaion  and  daclination  iam,  8a),  tha  ralationahip 
batwaan  tha  orbital  fraaa  and  tha  inartial  fraaa  (aaa  fig. 
3)  bacoaaas 


o1(t) 


coa(aa)  0  ain(aA)  01  0 


02(t)>»  0  coa(8fl)  -ain(8a) 


0  0  -1 


og(t)  0  ain(8a)  coa(8JI)  l-ain(aa)  0  coa(aA)  -10  0 


Multiplying  tha  aacond  and  third  aatricaa  yialdai 


r  ^  r  -i 

Oj(t)  10  0  -ain(aa)  coa(aa)  0 

og(t)  ►  *  0  coa(8a)  -ain(8a)  0  0-1 


o3(t) 


0  ain(8a)  coa(8fl)  -coa(aa)  -ai n(aA)  0  IK 


•J 


For  the  purpose  of  this  paper,  the  orbital  model  will 


consist  of  a  geocentric  satellite  in  two-body  motion.  In 


this  case,  8S  *  0,  and,  if  we  reference  all  time 


measurements  to  a  crossing  of  the  first  point  of  Aries 


(tg0),  then  we  can  solve  for  the  satellite's  right 


ascension  asi 


aa  «  fl(t-tgQ)  (24) 


It  is  noteworthy  that  the  attitude  prediction  to  be 


developed  herein  could  easily  handle  any  degree  of 


complexity  orbital  model,  so  long  as  the  satellite's 


epheaerides  were  known  quantities,  and  the  satellite’s 


orbit  produced  no  short-term  perturbations  which  night 


alter  the  orientation  dynamics  nodal.  In  the  simplified 


case  we  have  deve loped « 


°2  ?  m 


-slnC0(t-tgo) ]  cos(ft( t-tgQ) 3  0  I 


-coa[0<t-tgo>]  -sintn( t-tgo) ]  Ol  [K 


-1  J  (25) 


Rearranging,  the  inertial  coordinates,  in  terms  of  the 


orbital  frame,  becoaet 


I 


K 


-ain[fl(t-t 

coaCClCt-t 

0 


gO>  3 

0 

-cosCD(t-tg0) ] 

°1 

go*  ^ 

0 

-sin[0(t-tgo) ] 

< 

°2> 

-1 

0 

°3 
v.  J 

(26) 


Now,  given  a  known  star  position  (a^,8^>,  (aae  fig.  4) 
the  coordinates  can  be  transforaed  into  the  I  J  K  systems 


x,r 


Figure  4.  Star  Vector  Definition 


1^  =  coa(a^)coa(8^)  I+sin(a^ )cos(S^ )  J+sintS^)  K  (27) 


Now,  coabining  equations  (26)  and  (27)  with  the 
relation  between  the  orbital  axes  and  the  body  axes  yields: 


T 

•  " 

ll 

cos(S^ )cos(a^ ) 

-sin [D ( t-tg0 ) ]  0 

-cos[n(t-tgo) 3 

p 

bi 

l2 

>•< 

cos(8^ )sin(a^ ) 

► 

cosCD( t-ttfo ) ]  0 

-sinCn(t-ttfo)3 

R 

< 

b2 

x3 

sin(81) 

h  - 

0  -1 

0 

.b3 

(26) 
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where » 


[R]* 

COB ( tjlg  )  -Sin(cDg)  0 
sin((0g)  coa(tllg)  0 
0  0  1 
(28) 

Expanding  the  £R]  aatrlx  will  aake  it  nora  convenient 
for  use  later.  Tha  expansion  will  be: 

[R]  «  CRKaigjnRacaij^JlCRStOig)]  (30) 


cos((|>2)  0  sin((|>2) 
0  10 
-sin(t|>2)  0  cos(<J)2> 


10  0 
0  cos(4»i)  -sin(4>i) 
0  sin(Oi^)  cos(4>i) 


where < 


CRl(<p2)]* 


cos ( 4)2 )  0  sin(4>2> 

0  10 

I 

-s in(4)2>  0  cobI^q) 


(31) 


[R2(4>i )  3 


10  0 
0  cos(4)i)  -sin(4)j_) 
0  sin(4*i>  cos(4*x) 


(32) 


CR3(4>3)  3* 


cos(4>3)  -sin(4>3)  0 
sin ( 4>g  )  cos(4>3)  0 


(33) 


Star  Sanaor  Oparation 


Thia  papar  will  utiliza  a  v-alit  atar  aanaor  to  obtain 
atallar  data.  Although  a  major  portion  of  tha  oparation  of 
a  atallar  attituda  dataraination  achaaa  involvaa  atar 
idantif ication,  tha  aathod  to  ba  davalopad  harain  will 
aaauaa  a  priori  atar  idantif ication  with  aach  atar 
aighting.  Tha  atar  aanaor  will  ba  dafinad  to  hava  tha 
firat  alit  alignad  parallel  with  tha  bg  axi a,  with  tha  b^ 
axia  out  tha  boraaight  axi a,  and  tha  bg  axia  will  coaplata 
tha  daxtral  aat  (aaa  Figura  5).  Phyaical  aiaoriantationa 
could  ba  aaaily  accountad  for  (raf  7 >218-221),  but  thay 
will  ba  naglactad  in  thia  raport. 


Each  alit  of  tha  atar  aanaor  haa  a  planar  fiald  of 
vlaw.  Tha  aanaor  ia  fixad  to  tha  aatallita,  and  thua  acana 
tha  aky  am  tha  aatallita  apina.  Tha  width  of  tha  fiald  of 
viaw  will  ba  takan  am  thraa  dagraaa. 

At  a  tiaa  (t)  whan  a  atar  antara  tha  fiald  of  viaw  of 
alit  ona  and  ia  idantifiad,  tha  following  ralationa  axiat 
(aaa  Figura  6)s 


l^(t)  •  t>2  •  0  ("no”  arror)  (34a) 

li(t)  •  bj  «  coa(8li(t))  (34b) 

l*(t)  •  bg  ■  ain(811(t)>  (34c) 

$li(t)  -  1.5  dagraaa  (34d) 


Aftar  a  abort  tiaa  pariod  (tg),  tha  atar  will  noraally 
paaa  through  alit  two,  and  a  naw  gaoaatric  ralationahip 
axiata  (aaa  Figura  3). 

Thia  naw  ralationahip  can  ba  daacribad  in  taraa  of  a 
naw  body  axla  aat  (cj_,  c2,  eg)  amt 


lj(t>tg)  •  cq  ■  0  ("no”  arror) 

li(t+tg)  •  c^  »  coa(82i( t+tg) ) 

1 1 ( t+tg )  •  eg  ■  ain(82i ( t+tg ) ) 

®2i ( t+tg )  -  1.5/008(0!)  dagraaa 


(35a) 


(35b) 


(35c) 


n 


(36d ) 


Figura  6.  Star  Slit  Rafaranca  Fraaaa 

Tha  body  c^.cg.cg  fraaa  i«  linkad  to  tha  principal  body 
fraaa  through  two  anglaa,  0Z  and  6^,  which  ara  propart iaa 
of  tha  atar  aanaor’a  daaign  (aaa  Figura  7).  Froa  figura  7, 
tha  choica  of  tha  v-alita  in  tha  atar  aanaor  bacoaaa 
apparanti  tha  v-alita  allow  tha  daclination  of  tha  atar  in 
tha  body  fraaa  to  ba  datarainad.  Tha  ralationahip  batwaan 
cl»c2,c3  *n<*  b^,b2,bg 


coa(dz)  -ain(0z)  0  I  bj_ 


eg \z  0  com(Q^)  ain^)  ain(8z)  coa(8z)  0  b2>  (38) 


eg  0  -ain(e^)  coaCG^) 


1  bg 


which  raducaa  toi 


Figure  7.  Star  Scannar  Dasign  Anglaa 


f  ' 

C1 

cos ( ez ) 

-sin(0z) 

0 

bl 

<c2 

>« 

cos(61 )sin(8z) 

cos(0i )cos(0z) 

sin(01) 

H 

b2 

<=3 

-sin(0i)sin(0z) 

-sin(0i)cos(0z) 

cos(0i ) 
• 

b3 

Rafarring  to  aquation  (28),  latting  tha  star  position 
vactor  and  orbit  position  aatrix  ba  raprasantad  by  P(a^t 
8^  )  and  RQ(t-tgo),  raspactiva ly ,  than  coabining  this 
axprassion  with  aquations  (34), (36),  and  (37),  yialds: 


{p(a1,81)jr[Rn(t-tg0)][Rl(a»2)][R2(^i)][R3(M»3)]  fbj  »  Jpbli)  (38) 
{P(ai,8i)}T[Rn(t-tgo)][Rl(a»2)][R2(4»i)][R3(M>3>][R0]  *  {Pb2iJ  (39> 


wharat 


jpCaj.S^ 


co «(a^ )co«(Si  ) 
sin(a^ )cos(S^ ) ^ 
ain(8i ) 


(40) 


jRn<t-tgojj 


-sin[0(t-tgo)3  0  -coa[0(t-tgo) ] 
cos(fl(t-tgo)  3  0  -sin[Q(t-tgo) 3 
0-1  0 


(41) 


[••] 


COS(0Z) 

-sin(0z) 

0 

f 

bl 

cos(0^ )sin(6z) 

cos(9i )cos(8z) 

sin( ) 

< 

b2? 

-sin(0i )sin(0z) 

-sin(0i )cos( 0Z ) 

cos(0i ) 

b3 

b  - 

(42) 


coa(8li( t ) ) 
0 

ain(Sli(t)) 


(43) 


Icoa(82i( t+ts) ) 
0 

flin(S2i< t+ts) ) 


(44) 


A a  a  final  siaplif ication,  let  the  left  hand  side  of 
equations  (38)  and  (39)  be  represented  by  Obsl  and  0bs2, 
respectively.  Then,  using  only  the  equality  portions  of 


these  equations,  the  observation  relations  reduce  to> 


0bs22  *  0 


(45b) 


Truth  Model 

An  operating  star  scanner  would  obtain  inputs  from 
stellar  slit  transits,  which,  after  being  tiae-tagged  by 
the  satellite’s  reference  clock,  would  be  sent  to  the 
attitide  deteraination  package.  These  slit  crossing  tiaes 
would  be  directly  dependent  on  the  satellite’s  dynaaic 
characteristics.  In  computer  simulations,  the  observation 
inputs  must  be  provided  by  the  truth  aodel.  The  model  used 
in  this  effort,  as  stated  earlier,  assumes  a  satellite  in 
geosynchronous,  two  body  orbit,  with  nominally  earth¬ 
pointing  attitude.  Under  these  assumptions,  the  reference 
attitude,  (the  orbital  reference  frame),  is  established  by 
the  time  elapsed  since  tg0,  and  this  reference  rotates 
about  the  o2  axis  at  constant  rate  (0)  of  7. 292115858xl0-5 
rad/sec  (10:429).  Thus,  given  the  true  initial  attitude 
state  parameters  (w10,  w2q,  W30*  ^10*  ^20*  ^30 *  and  the 
initiml  time,  with  respect  to  tgQ,  the  true  state  of  the 
satellite  can  be  tracked  by  numerically  integrating  the 
exact  first  order  differential  equations.  Knowledge  of  the 
true  state  permits  the  star  vector  to  be  transformed  into 
the  star  sensor  slit  frame,  using  the  transf ormations 


developed  in  the  last  section.  Watching  the  b2<C2) 
component  of  a  star  vector  for  a  sign  change  allows 
determination  of  a  slit  plane  crossing  time.  After 
checking  f ield-of-view  constraints,  a  star  sighting  can  be 
recorded. 

For  this  research,  a  star  data  base,  which  is  accessed 
by  the  truth  model  algorithm,  was  compiled  from  the  Bright 
Stmrs,  J1984.6  table  of  "The  Astronautica 1  Almanac”  (4iHl- 
H31).  This  data  base  lists  a  circumpolar  band  of  stars, 
approximately  30°  wide,  centered  on  a  nominal  satellite 
position  of  0°  right  ascension  at  tgo*0  (see  figure  8). 

The  list  contains  301  stmrs  of  visual  magnitude  +6.0  or 
brighter  (see  Table  1,  Appendix  A). 


XU 


(  Sat.  boresight 
tor  tQ0  =0 ) 


Figure  8.  Truth  Model  Star  Field 


Estimator 


Tha  final  component  of  the  attitude  determination 
package  is  the  estimator.  Here,  the  truth  data  is  combined 
with  the  approximated  dynaaios  and  observation  relations  to 
fora  the  state  estimate.  Due  to  the  nonlinear  nature  of 
the  dynamics  and  observations  relations,  some  nonlinear 
estimation  algorithm  is  desired.  The  choice  made  for  this 
research  is  a  nonlinear  least  squares  algorithm. 

Nonlinear  Least  Squares  Algorithm 

A  typical  nonlinear  least  squares  algorithm  operates  in 
the  following  manners 

I)  Propagate  the  state  vector  to  the  observation 
time.  Also,  obtain  the  state  transition  matrix  *(t,t0). 

II)  For  each  observation,  calculates 


a)  rA  *  z*  -  G(xr#£( ti ) , ti ) 


b)  H±  *  OG/ax)  C(xre£(t1), 


c)  *(titto)  a  dx< tA )/8x(t0) 


d)  Ti  =  Hi*1(ti,t0) 


(46) 


(47) 


(48) 


(49) 


III)  Add  new  terms  to  the  running  sums  of  the  aatrixs 


and  the  vector: 


IV)  Compute  the  covariance  of  the  correction: 


P8x  *  CS  (60) 


and  the  state  correction  at  epoch: 


8x(t0)  *  P8x  S  TiTQ^ri  (61) 


V)  Correct  the  reference  attitude  state: 


xre£(new)  *  xre£ (old)  +  Sx(tQ)  (62) 


VI)  Repeat  steps  I  through  V  until  convergence  is 
achieved.  (Check  residuals  for  valid  convergence.) 

(11:68  &  12:60,31) 

The  components  needed  to  employ  this  algorithm  will  now 


be  described 


The  observation  relations  developed  earlier  <eq.  45), 


when  expressed  as  a  vector,  make  up  the  observation  vector 
<6(x^,t^>,  which  can  be  expressed  ast 


!Ob«l(x^,  tn  >  1 

r 

0bs2(xi2» t12)J 

(63) 

Some  simplification  can  be  obtained  by  treating  the 
information  from  each  star  scanner  slit  as  a  separate  data 
point.  In  this  manner,  the  data  vector  {z>  is  reduced  to  a 
scalar  value,  and  the  observation  vector  CG)  relation  is 
reduced  to  two  scalar  relations.  Which  relation  is  used  is 
determined  by  the  slit  making  the  observation. 


Gradient  of  the  Observation  Vector  [H] 

Now,  having  tranforaed  the  observation  relation  into  a 
scalar  relation,  the  gradient  of  the  observation  relation 
with  respect  to  the  state  vector,  the  CH]  matrix,  becomes  a 
row  vector  relation.  In  order  to  obtain  CH],  the  following 
partial  derivatives  will  be  required! 

(64) 


d/dtp2  [Rl(<]>2)] 


-sin(«i>2)  0  cos((i>2) 

0  0  0 
-cos(il>2)  0  -ain(<|J2) 


d/di^  CR2((Ji^)3  «  0  -*in(4»i)  -co «<4»i)  (65) 

0  co«((^i)  -cin(<t>^) 

-•in ( (tig)  -co*((Hg)  0 

6/ditig  C R3 ( (Jig )  ]  *  co* ((Jig)  -*in(<pg)  0  (56) 

0  0  0 

With  raapact  to  tha  atata  variablaa,  *11  othar  partial 
darivativaa  art  ztro,  Aaaigning  aatrix  noaanclatura  to  tha 
abova  aquation*,  lati 

[HR1((|>2)]  *  6/6(Jig  tRl<M»2>3  (57) 

CHR2(iH1)]  *  6/6*3,  CR2((Hi)]  (56) 

CHR3(0ig)]  »  6/6*g  [R3(*g)3  (60) 

Now  CH3  can  ba  aaaaablad  ast 

CH3  *  [  hg  hg  h4  hg  hg  3  (60) 

Tha  alaaant*  of  (H3  arat 


h^  *  dS/dwj_ 

(61a) 

hg  *  66/ dug 

(61b) 

hg  *  66 /dug 

(61c) 

hg  *  dG/d«t>2 
hg  *  dG/3tp3 


(61e ) 


(61f ) 

Sine*  the  observation  relations  ere  dependent  on  the 
observation  slit,  the  eleaents  of  [H]  will  also  vary  with 
the  observation  slit.  For  slit  1  star  sightings,  the 
eleaents  of  CH]  aret 


hj  *  0  (62a) 

h2  ■  0  (62b) 

h3  *  0  (62c) 

h4  ■  El2  of  {{P(a1,6i)>TCR0(t-tgo)nRl(4J2)] 

*  CHR2(a»1)]CR3(M>3)]>T  (62d) 

h6  *  El2  of  ((P(ai,6i))TCRn(t-tffo)nHRl(M)2)] 

*  CR2(<H1)][R3(0i3)]}t  (62e) 

he  *  El2  of  ((P(ai,Si)>TCRn(t-tgo)  ]  CRl(a>2)] 

*  CR2(Hi1)]CHR3(a>3)]>T  (62f ) 


For  slit  2  star  sightings,  the  eleaents  of  [H]  are: 


hj  *  0  (63a) 

(63b) 


h2  «  0 
h3  ■  0 

h4  »  El2  of  (CP(ai,Si))TCRn(t-t-0)HHRl(a)2)] 


(63c) 


(83d) 


*  cr2(«j>1)] [R3c«gmR<eit ez>]>T 
h6  »  El2  of  <<P(ai,8i)>TCRn(t-fctfo)][Rl(«|»2)3 

*  [HRR(*1)UR3(*g)HR<ei,ezH}T  (83e) 

he  «  El2  Of  ({Ptai.SiJjTCRfltt-tgo)]  CR1(0>2)] 

*  CR2(M>1)]CHR3(a>3)3[R(ei,0z)])T  (83f) 

Proa  equations  (80-63),  it  can  be  seen  that  the 
elements  of  [HI  for  slit  1  observations  are  determined  from 
element  2  of  a  vector,  which,  if  that  vector  is  post- 
multiplied  by  the  CR0]  matrix,  yields  the  elements  of  CH] 
for  slit  2  observations.  Problems  which  sight  arise  from 
the  small  difference  in  observation  times  of  a  particular 
star  by  each  star  sensor  slit  are  thus  avoided  by  treating 
each  sighting  as  separate  data.  All  that  is  now  required 
to  form  the  residual  is  the  transf ormation  of  the 
observation  relation  into  a  relation  which  estimates  a  slit 
crossing  time. 

Computing  the  residua  1  (r) 

The  truth  model  gives  the  time  for  a  given  slit 
crossing.  This  parameter,  tact,  permits  solution  of  the 
state  estimate.  In  order  to  compute  the  residual  (tacfc- 
t#at),  it  is  necessary  to  determine  the  estimate  of  the 
slit  crossing  time  from  the  estimated  state.  Since  t  is  an 


implicit  variable  in  the  dynamics  mnd  in  the  observation 


relation,  soae  iterative  acheae  will  be  required  to  find 
test*  The  observation  relatione  lend  theaeelvee  readily  to 
a  Newton-Rhapheon  iteration  technique,  and  this  acheae  will 
now  be  developed.  The  fora  of  each  observation  relation 
iet 

tf^(x,t)  *  0  (i«l,2)  (84) 

The  Newton-Rhaphson  aethod  to  solve  for  t  will  bet 

t„ew  *  told  “  tfi/<<S*i/dt>  (85) 

The  actual  slit  crossing  tiae  will  be  used  as  the 
initial  seed  for  each  iteration.  In  order  to  find  the  tine 
derivatives  of  each  observation  relation,  the  chain  rule 
will  be  applied,  as  follows: 

dgA(x,  t  )/dt  *  dgj/dt  +  (dgj/dxHdx/dt)  (88) 

The  only  new  expressions  needed  here  are  the  partial 
derivatives  of  the  observation  relations  with  respect  to 
tiae.  Recalling  the  observation  relations  (eqs.  38,88)  the 
only  coaponents  of  these  equations  which  are  explicit  in  t 
are  the  eleaents  oft 


[Rn<t-fcffon 


The  partial  darivativa  of  this  axprassion  with  raspact 


to  t  ist 


cos[fl( t-tgo ) ]  0  -sinCfl(t-tgo) ] 


3/atCRn(t-tgo)]  =  -0  sin[fl(t-tgo)]  0  cos[0(t-tgo)]  (67) 


Now,  tha  first  taras  on  tha  right-hand  sida  of  aquation 
(66)  becoaai 


Elg  of  {(P(a1,61)>TCd/dt(Rfi(t-tgo)]][Rl(<|>2)] 

*  [R2(4)1)][R3(  (Jig  )  3  >  (68) 

El2  of  UP(a1,S1)>TCd/dt[RO(t-tgo)]]  [R1(0»2>1 


*  [R2(O)1)HR3(Mi3)HR(0i,ez)]} 


(69) 


The  other  two  ter as  required  are  (dg/dx)  and  (dx/dt), 
which  are  already  developed  as  CH]  and  the  first  order 
state  differential  aquations.  Now,  aquation  (66)  becomes t 


dg/dt^  *  g^  +  CH^J  (dxj/dt^) 


(70) 


where  i(*l,2)  signifies  tha  applicable  star  sensor 


illt. 


WWW 
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Subatitution  of  aquation  <70)  into  aquation  (66)  yialda 


tha  baaic  aquation  of  tha  Nawton-Rhaphaon  achaaa,  which, 
aftar  itarating  to  a  apacifiad  tolaranca,  givaa  an 
aatiaation  of  tha  alit  croaaing  tiaa.  Tha  raaidual  acalar, 
r^,  can  ba  calculatad  aat 

►  ri  *  fcact*.  “  taat*.  <71) 

Data  Sradiant  CT] 

f  Having  found  tha  raaiduala,  tha  gradiant  of  tha  data 

acalar  with  raapact  to  tha  initial  atata,  CT],  ia  tha  only 
raaaining  ralation  raquirad  for  tha  aaaambly  of  tha 
*  aatiaata  corraction  covarianca.  Tha  darivation  of  CT]  ia 

aa  followat 


ri  a  ( fcac t*. -feast *•>  *  tHi]  C8xiCti)>  (72) 

wharat 


CBXiCti)}  *  C*(tift0)]  C8x(t0)>  (73) 


Subatitution  of  aquation  (73)  into  aquation  (72) 
yialdai 


rA  «  CHi)  C*(titt0)]  {8x( t0 ) )  (74) 


The  C T ]  matrix,  (herein  a  row  vector),  is  simply  a 
shorthand  notation  for  the  product  of  C H ]  and  t*]i 

[Til  *  CHi]  C*(t1,t0)]  (76) 

In  terns  of  the  data  scalar  and  the  state,  the  [T] 
matrix  can  be  considered  to  bet 

*  atj/dxQ  (78) 

The  problem  of  t  being  an  implicit  variable  in  the 
observation  relations  now  resurfaces.  This  problem  is 
circumvented  by  obtaining  the  [TJ  matrix  numerically  using 
a  finite  differencing  technique.  For  each  observation  tine 
t^,  an  estimated  slit  crossing  time  is  computed.  By 
varying  each  initial  estimated  state  variable  by  a  small 
amount,  a  new  estimated  slit  crossing  time  can  be  obtained 
for  each  variable  change.  The  approximate  relationship 
which  results  isi 

[T±]  X  Atj/Axo  (77) 

Use  of  the  approximate  relationship  for  [T]  eliminates 
the  need  to  compute  the  [#]  matrix. 


Data  Covariance  Matrix  CQ3 

The  data  covariance  matrix  is  a  property  of  the 
satellite  hardware.  It  is  a  numerical  declaration  of  the 
precision  of  the  instrument's  ability  to  record  the  data. 

A  simplified  approach  is  taken  in  this  problem  by  treating 
the  data  covariance  as  a  scalar  value.  Thus,  the  value  of 
q  will  simply  become  the  square  of  the  standard  deviation 
of  the  observation  time  for  a  particular  start 

^i  *  <at*.)2  (78) 

Now,  all  components  required  to  compute  the  covariance 
of  the  estimate  correction  have  been  derived. 

Covariance  o£_  the  estimate  correction  CP ] 

The  final  component  of  the  least  squares  estimation 

algorithm  to  be  computed  is  the  covariance  of  the  estimate 
correction,  or,  as  it  is  commonly  called,  the  [P]  matrix. 
As  shown  by  equation  (50),  the  [P]  matrix  is  obtained  by 
inverting  the  matrix  described  by: 

CTiTQi-lT1] 

With  the  aodif ications  described  above,  the  resulting 
expression  for  the  C P 3  aatrix  in  this  simulation  is: 


VI)  Vary  each  initial  state  estiaate  variable,  and 


VII) 


VIII) 


IX) 


repeat  step*  III-V  to  obtains 

T* 'Z.  A ti/Ax0 

Add  new  tarae  to  the  running  suae  of  the  aatrixs 

q-lCST^Til 

and  the  vectors 

q^CSTi1^] 

Coaputa  the  covariance  of  the  corrections 

P8x  «  q  [  [ST^Ti  ]”1 3 

and  the  state  correction  at  epochs 

8x(tQ)  »  qCPsxHSTj^ri) 

Correct  the  estiaated  attitude  states 

xtat(new)  *  xttflt(old)  ♦  8x(t0) 


,n 


Chapter  Three  Program  lap leaentation 


A  major  portion  of  the  effort  in  thia  thesia  involvea 
the  creation  of  a  working  aiaulation  of  a  aatellite 
attitude  atate  prediction  aodel.  Thia  chapter  ia  devoted 
to  taking  a  brief  qualitative  look  at  the  prograa a  uaed  in 
the  aiaulation.  Alao  diacuaaed  are  aoae  of  the  hardwired 
conatrainta  uaed  in  the  prograa,  and  the  reaaoning  behind 
the  conatraint  choicea  (where  pertinent).  For  a  coaplete 
prograa  liating,  aee  Appendix  B. 


Prograa  Overview 

The  prograa  uaed  to  aiaulate  the  aatellite  attitude 
eatiaation  conaiata  of  a  aain  controller  prograa,  16 
aubroutinea,  and  two  input  filea.  All  prograa  code  ia 
written  in  FORTRAN  77  language,  and  one  IMSL  aubroutine 
(LEQT2F)  (ref  13tLEQT2F),  a a  well  a a  one  IMSL  function 
(G6NQF)  (ref  13>GGNQF) ,  are  utilized. 

Input  Fi lea 

There  are  two  input  filea  for  the  eatiaation  prograa,  a 
Bright  Star  Catalog,  "atardat”  (App  B2&-B30),  and  a  initial 
conditiona  input  file,  ’’gueaain”  (App  B14).  The  Bright 


Star  Catalog  containa  301  atara,  liating  each  atar’a  right 
aacenaion  and  declination  (in  radiana),  a a  well  a a  the 


star's  identification  nuaber  and  visual  magnitude 


and  a  noisa  aaad  for  tha  IMSL  function’*  paaudo-randoa 


noraal  nuabar  generator. 

Main  Program  Operation 

Tha  controlling  program  for  tha  attituda  aatiaation 
package,  ’’estimator,  f "  (saa  app  B2-B5),  reads  tha  initial 
conditions  input  file,  controls  tha  overall  program  flow, 
and  produces  tha  desired  output  information  in  a  usable 
format.  After  reading  tha  initial  conditions  input  file, 
tha  initial  conditions  are  stored  in  array  format,  and  are 
printed  to  tha  output  file.  For  each  epoch  update,  tha 
following  processing  occurst 

1)  Input  data  for  an  estimation  sequence  are 
initialized. 

2)  Prograa  control  is  transferred  to  subroutine 
"truth”  (App  B36-B37),  where  the  truth  aodel  star 
sighting  tiaes  for  on*  scan  period  are  determined  froa 
true  initial  conditions,  stellar  data,  observation 
relations,  and  propagation  of  the  true  state.  Star 
sighting  information  is  stored  in  a  temporary  data  file 
"sitings”  (App  01-D2),  and  program  control  is  passed 
back  to  the  aain  prograa. 

3)  Prograa  control  is  passed  to  subroutine  ’’torder" 


(App  B32-B33),  which  reads  the  ’’sitings"  file,  adds 


noise  to  the  truth  aodel  data' (if  required),  prints 
noise  inforaation  to  tha  output  file,  tiae-orders  tha 
star  sighting  inforaation,  and  storas  tha  adjusted 
sighting  inforaation  in  taaporary  data  fila  "stars" 
(App  D3-D4).  Prograa  control  is  than  passed  back  to 
tha  aain  prograa. 

4)  Next,  prograa  control  is  passed  to  subroutine 
"tresid"  (App  B34),  where  a  Newton-Rhaphson  iteration 
procedure  is  used  to  coapute  tha  estiaated  star 
observation  tiaes.  True  star  sighting  tiaas, 
observation  relations,  approxiaate  dynaaics,  and  first 
order  exact  state  differential  equations  are  used  to 
obtain  results.  True  observation  tiaes,  estiaated 
observation  tiaes,  and  stellar  inforaation  are  stored 
in  taaporary  data  file  "statfile”  (App  D6-D7).  Progra 
control  is  then  passed  back  to  the  aain  routine. 

5)  Prograa  control  is  then  passed  to  subroutine 
"dtdxO"  (App  BIO),  where  a  finite  differencing 
technique  is  used  to  coapute  the  CT1  row  vector  for 
each  star  sighting.  The  procedures  explained  in  the 
previous  step  are  used  to  obtain  results.  Star 
inforaation,  true  and  estiaated  slit  crossing  tiaes, 
and  CT]  data  are  stored  in  teaporary  data  file  "dtfile 
(App  08-09),  and  control  is  passed  back  to  the  aain 


prograa 


Finally,  after  tha  prograa  haa  coaplafcad  all  estiaation 
propagation  cyclaa,  aubroutina  "plotarr"  (App  B24)  raada 
tha  arror  data  fila  "err list",  puta  thia  inforaation  into 
plottabla  foraat,  and  atoraa  tha  inforaation  in  data  fila 
"errplot”  (D14-017). 

Othar  Routinaa 

Bafora  discussing  tha  hardwired  paraaatara  uaad  in  thia 
aiaulation,  tha  aubroutinaa  which  hava  not  already  bean 
specifically  discussed  will  be  presented. 

Subroutine  "theslrhs”  (App  B31)  contains  tha  exact 
first-order  differential  state  aquations.  Thia  routine  is 
accessed  by  tha  nuaarical  integration  aubroutina,  and  by 
tha  aubroutina  which  processes  tha  Wewton-Rhaphson 
iteration  for  the  aatiaatad  star  observation  tiaa. 

Subroutine  "haaing"  (App  B15-B17)  is  a  fourth-order 
predictor-corrector  nuaarical  integration  routine.  It  is 
uaad  to  propagate  tha  true  state  aquations  forward  in  tiaa. 
Thia  algoritha  waa  developed  by  Professor  Anderson,  Harvard 
University,  in  the  1980 'a.  The  version  used  in  this 
aiaulation  was  provided  by  Dr.  Williaa  Wiesel,  Professor  of 
Astronautics  1  Sciences,  AFIT. 

Subroutine  ”dyno2"  (App  B12,B13)  contains  the  second 
order  approxiaate  solution  to  the  state  attitude  dynaaics. 
Given  an  initial  tiaa,  state,  and  tiaa  of  concern,  this 
routine  coaputes  the  approxiaate  state  at  the  tiae  of 


concern.  A  aiailar  program  "dynol"  (App  Bll),  contain#  the 
firat  order  approximate  eolutione  to  the  etate  attitude 
dynamic#.  With  very  ainor  program  aodif icatione,  the 
eetimation  routine  can  be  run  with  firet  order  approximate 
dynaaica  in  place  of  the  current  aecond  order  model. 

Subroutine  "obaerv2"  (App  B20-B23)  contain#  the  matrix 
relation#  required  to  coapute  the  obaervation  relation# 
(g±)t  the  CHI  row  vector,  or  the  partial  derivative  of  the 
obaervation  relation  with  reapect  to  time  (gdot^).  The 
calling  routine  provide#  atar  information,  tine,  aidereal 
tine,  atate,  and  a  apecif ication  flag.  The  apecif ication 
flag  aignala  thia  routine  which  relation  ia  required. 
wobeerv2M  utilize#  two  matrix  operation  routine# i 
"aatxaat”  (App  BIB),  which  aultipliea  two  aatricea,  and 
"vecxnat”  (App  B40),  which  preaultipliea  a  matrix  by  a  row 
vector.  The  matrix  manipulation  routine#  were  provided  by 
fellow  AFIT  atudent  Capt  Keith  Greer,  GAE/88D. 

Subroutine  "tatate"  (App  B38,B39)  drive#  the  Newton- 
Rhaphaon  iteration  acheae  which  ia  uaed  to  compute  the  CT] 
row  vector. 

Subroutine  ”noiae"  (App  B19)  uaea  IMSL  function  ggnqf 
to  provide  a  one-aigma  paeudo  random  normal  time  error. 

Thia  time  error  ia  added  to  the  true  alit-croaaing  tine  to 
generate  "noiay”  truth  data  where  required. 


Hardwired  Parautetri 

Several  parameter  constant a  art  aabaddad  within  tha 
•atiaator  programs,  but  could  ba  altarad  with  ainor  program 
aditing.  Thaaa  altarabla  conatanta  will  now  ba  diacuaaad. 

Tha  aiza  of  tha  input  atar  data  baaa  ia  fixad  at  301 
atara  apanning  alaoat  a  30°  wida  longitudinal  awath  of  tha 
aky.  Incraaaing  tha  aiza  of  thia  liat  would  raquira 
adition  of  tha  STAR. OB  tabla,  but  tha  liat  could  ba  aaaily 
diainiahad  by  aditing  aithar  tha  aain  program  or  "atarcnv. 
f”  to  diacriainata  according  to  viaual  aagnituda. 

Tha  aatiaation  paokaga  worka  with  an  undar lying 
aaauaptlon  of  a  noalnal  aatallita  apin  rata  of  B 
ravolutiona  par  alnuta.  A  single  aatiaation  loop  procaaaaa 
all  data  obaarvad  in  a  typical  aingla  acan  period  (12 
aaconda).  In  ordar  to  raduca  tha  procaaaing  tiaa  involved 
in  generating  tha  truth  aodal  data,  tha  truth  aubroutine 
waa  prograaaad  only  to  look  for  atar  aightinga  for  a  tiaa 
period  equal  to  currant  epoch  tiaa  +  12  aaconda.  Also, 
after  a  atar  observation  is  recorded,  no  second  sighting 
for  that  atar  ia  permitted.  Altering  tha  typical  spin  rate 
of  tha  aatallita  would  raquira  changing  tha  parameter 
"tend"  in  tha  aain  program  to  a  value  approximately  equal 
to  tha  typical  acan  period  (in  aaconda).  Altering  tha 
truth  aodal  operation  further  would  raquira  sore  extensive 


program  change*,  and  this  topic  will  ba  discussed  latar  in 
this  report. 

Star  scanner  physical  parameter*  0^  and  0Z  are 
currently  set  to  15°  and  0.5°,  respectively,  and  they 
reside  in  subroutines  ”ob*erv2”  and  "truth”.  The 
characteristics  of  the  star  scanner  could  be  altered  by 
siaply  editing  these  paraaeters.  The  scanner’s  field  of 
view,  also  hardwired,  is  set  to  three  degrees  in  subroutine 
"truth”.  Changing  the  parameter  ”fov”  in  "truth”  is  all 
that  is  required  to  alter  the  scanner’s  field  of  view. 

The  orbital  rate  paraaeter,  0,  carries  with  it  many 
implications.  In  essence,  this  parameter  suaamrizes  the 
assumptions  of  a  satellite  reference  attitude  aodel  which 
arises  froa  a  circular,  two-body  orbit  having  0° 
Inclination.  It  also  might  appear  to  iaply  that  the 
reference  attitude  maintains  its  nominal  earth-pointing 
state  without  torques.  This  is  not  the  case,  however,  for 
although  the  governing  state  equations  assuae  torque-free 
motion,  the  estimation  algorithm  could  easily  handle  the 
presence  of  small,  impulsive,  attitude  correction  torques 
(so  long  as  the  true  state  corrections  were  presented  to 
the  truth  aodel).  Analysis  of  a  more  complex  equatorial 
reference  attitude  aodel  would  require  replacing  the  0 
paraaeter  in  routines  "theslrhs"  and  ”observ2”  with  a 
subroutine  which  described  the  reference  attitude  rotation 


rata.  Tha  study  of  inclinad  orbits  would  require  a  revised 
derivation  of  the  observation  relations. 

Residual  rejection  criteria  are  established  in 
subroutine  "covarian".  Currently,  these  paraaeters  are  set 
to  3  seconds  for  the  first  two  estiaation  iterations,  and  3 
sigaa  for  succesive  iterations.  If  all  residuals  are 
rejected  on  one  pass,  the  constraints  are  widened  by  an 
order  of  aagnitude.  For  different  case  studies,  it  sight 
be  desirable  to  edit  this  rejection  scheae. 

The  finite  differencing  technique  used  in  this  effort 
uses  a  forward-differencing  approach,  with  the  initial 
state  being  varied  by  a  finite  difference  aaount  of  0.1%. 
Probleas  aight  arise  for  cases  where  initial  conditions 
approach  zero.  Different  criteria  could  easily  be  arranged 
by  saall  editions  to  the  "covarian”  subroutine. 

For  this  research,  the  inertia  ratio  ((A-C)/A)  is  set 
to  a  value  of  -1.  This  represents  a  typical,  stable,  tuna- 
can  type  of  satellite  with  a  spin  aoaent  of  inertia  equal 
in  aagnitude  to  twice  the  aagnitude  of  the  transverse 
aoaents  of  inertia.  To  study  a  different  case,  the 
paraaeter  "ak"  in  the  "thelrhs”  and  "dyno”  routines  would 
have  to  be  edited. 

In  the  subroutine  "truth",  a  siaplif ication  is 
incorporated  which  locks  sightings  Bade  by  slit-2  of  the 
star  scanner  to  slit  1  observations.  Essentially  the 


f ield-of-view  of  a lit  2  is  adjuated  ao  that  all  atars 
obaervad  by  a lit  1  ara  alao  obaarvad  by  a lit  2.  For  a  caaa 
atudy  involving  larga  aatallita  pracaaaion  rat aa,  thia 
critarion  night  naad  to  ba  adjuated. 

Tha  final  hardwirad  conatrainta  uaad  in  this  aiaulation 
ara  tha  tolarancaa  in  tha  nuaarical  intagration  routina 
"haiing",  and  tha  Newton-Rhaphaon  convarganca  (aubroutine 
"tatata").  Whila  thaaa  tolarancaa  provad  adequate  for  thia 
caaa  atudy,  it  night  ba  daairabla  to  verify  thaaa 
paranatara  in  a  different  analyaia. 

Now  that  a  full  deacription  of  tha  progran  operation 
and  theoretical  aaaunptiona  have  bean  preaented,  raaulta 
froa  aavaral  caaa  atudiaa  can  ba  preaented. 


Chapter  Four  Results 


This  section  includes  results  froa  3  case  studies,  as 
well  as  results  obtained  while  verifying  some  of  the 
subroutine  algorithms.  During  this  presentation, 
qualitative,  as  well  as  quantitative,  results  are  compiled. 
The  sequence  used  to  present  the  results  parallels  the 
sequence  utilized  to  develop  the  siaulation  program. 

Typica 1  Input  Fi lea 

For  the  case  studies  which  will  follow,  two  typical 
cases  were  considered.  For  the  first  case,  all  true 
initial  conditions  were  set  to  zero  except  for  the  spin 
rate  (5  revolutions  per  minute),  and  the  initial  spin  angle 
(0.8  radians).  This  case  was  chosen  to  depict  a 
precession-free  situation  which  should  produce  predictable 
results.  For  the  second  case,  it  was  desired  to  find  a 
general  Initial  true  state  which  contained  arbitrary 
(small)  deviations  froa  the  first  case.  To  achieve  this 
goal,  the  following  true  initial  conditions  were  chosen: 


*10 

z 

0.01 

rad/sec 

(0.67293° /sec) 

x20 

z 

0.06 

rad/sec 

(2. 88479°/sec ) 

x30 

z 

0.6223608773 

rad/sec 

( 29. 9290°/sec ) 

x40 

z 

0.06 

rad 

(2.88479°) 

X50 

z 

0.06 

rad 

(2.88479°) 

xeo 

z 

0.8 

rad 

(46.8366°) 

Eat  lasted  initial  condition*  ware  chosen  to  be  close  or 


equal  to  true  initial  conditions. 

Dynaaic*  Checkout  Package 

The  first  coaponent  of  the  attitude  estiaation  package 
to  be  verified  was  the  dynaaic*  algorithm.  The  results  for 
a  single  scan  run  with  case  two  true  initial  conditions 
show  that,  qualitatively,  both  the  first  order  and  second 
order  approximate  dynaaic  models  closely  matched  the 
numerically  integrated  exact  equations.  Graphs  of  the 
d>2«  and  <pg  states  are  shown  in  Figures  10,  11  and  12. 
(Results  from  w  states  are  not  shown  since  these  equations 
are  solved  exactly  in  closed  fora). 

While  the  single~scan  dynamics  plots  showed  close 
agreement  between  approximate^ and  truth  model  dynamics 
solutions,  they  did  not  clearly  show  the  second  order 
approximate  equations  presented  any  advantage  over  the 
first  order  approximate  relations.  For  that  reason,  the  * 
relations  were  next  plotted  in  a  similar  fashion  for 
Identical  initial  conditions.  Examples  of  these  results, 
shown  in  Figures  13,  14  and  15,  clearly  show  the 
superiority  of  the  second  order  approximate  dynamics  over 
the  first  order  approximations. 
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Taat  Ca an  Raaulta 

Of  tha  night  taat  c aaaa  which  warn  analyzad,  tha  firat 
four  caaaa  containad  no  noiaa  in  tha  obaarvation  data.  For 
tha  raaainlng  taat  caaaa,  paaudo-randoa  normal  noiaa  waa 
addad  to  tha  truth  nodal  alit  croaaing  tiaaa.  Tha  inducad 
noiaa  haa  a  aean  of  0  and  ona  atandard  daviation  aqunl  to 
3.1622S  *  10”®  aaconda  (  1  aigaa  ). 

Taat  Caaa  1  waa  run  to  analyza  a  aingla  acan  aatiaation 
with  caaa  1  trua  initial  conditiona.  Tha  initial  atata 
aatiaata  waa  aat  to  tha  following  valuaat 


9 

0.011 

rad/aac 

9 

0.066 

rad/aac 

w30  <  aat ) 

9 

0.6723698778 

rad/aac 

*10(nat) 

9 

0.066 

rad 

*20(*-t) 

9 

0.066 

rad 

*30 (aat) 

9 

0.808 

rad 

Tha  aatiaation  packaga  convargad  aftar  four  itaration 
cyclaa.  Tha  convargad  aatiaata  of  tha  initial  atata 
raducad  tha  total  pointing  arror  from  4.48  dagraaa  to 
10.336  arc  aaconda.  Total  pointing  arror  ia  dafinad  by  tha 
aquation t 

TPE  *  C(6*^q)^  C84>20^  *  C70) 


A  qualitativa  rapraaantation  of  tha  convargad  aatiaata 


is  shown  in  Figure  16.  This  graph  dspicts  a  on#  scan 
propagation  of  tha  trua  and  astiaatad  dynamics  of  tha  star 
scannar  borasight  axis  against  tha  background  star  fiald. 
Thasa  rasults  wara  fait  to  ba  quita  adaquata  for  a  fina 
attituda  astiaator. 

Tast  Casa  2  is  a  rapaat  of  tast  casa  1,  but  with  casa  2 
trua  initial  conditions.  This  casa  convargad  to  a  total 
pointing  arror  of  6.736  aro  saconds.  Rasults  from  this 
casa  ara  plottad  in  figura  17. 

Tast  Casa  3  was  tha  final  noise-free  sing la  scan 
astiaation  run  parforaad.  For  this  casa,  conditions  wara 
rapaatad  froa  a  starting  valua  of  0  to  -50  ainutas,  thus 
raworking  Tast  Casa  2  against  a  diffarant  star  fiald.  Tha 
convargad  total  pointing  arror  achiavad  was  40.08  arc 
saconds.  Rasults  froa  this  casa  ara  plottad  in  Figura  18. 

Tast  Casa  4  chacks  tha  aultipla  scan  capability  of  the 
astiaation  package.  For  this  casa,  initial  conditions 
identical  to  Tast  Casa  1  wara  used.  After  achieving  a 
convargad  state  estimate  at  epoch,  tha  approximate  dynamics 
wara  than  propagated  for  approximately  30  scans  (5  ainutas) 
before  another  astiaation  sequence  was  parforaad.  This 
process  was  rapaatad  for  six  astiaation  sequences. 

Rasults,  showing  total  pointing  arror,  ara  plottad  in 
Figura  18.  In  this  figura,  it  can  ba  seen  that  tha 
convargad  astiaata  for  tha  fifth  astiaation  sequence  is 
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aubatantia 1 ly  leas  accurate  than  the  raaulta  achiavad  with 
tha  other  aatiaatlon  aaquancaa.  Thia  reduced  accuracy  ia 
due  to  the  fact  that  only  5  atara  were  obaerved  during  the 
acan  uaed  to  generate  thia  eatiaation  aequence,  whereaa  the 
other  aequencea  proceaaed  at  leaat  10  atara. 

The  remaining  teat  caaea  depicted  realiatic  aiaulationa 
with  noiae  preaent  in  the  obaervation  data.  Teat  Caae  6 
repeata  Teat  Caae  1,  but  with  one  eigna  random  normal  noiae 
added  to  the  true  obaervation  timea.  Total  pointing  error 
for  thia  caae  converged  to  66.14  arc  aeconda  (aee  Figure 
20).  Teat  Caae  6  ia  a  repetition  of  Teat  Caae  1  with  noiae 
added.  The  total  pointing  error  for  thia  caae  converged  to 
67.16  arc  aeconda  (aee  Figure  21).  In  thia  figure,  the 
propagation  of  the  initial  eatimate,  aa  well  a a  the 
converged  eatimate,  are  diaplayed. 

Teat  Caae  7  inveatigatea  the  multiple  acan  capability 
of  the  eatiaator  with  ahort  propagation  perioda  (noiae 
preaent).  For  thia  caae,  10  eatiaation  eye lea  were 
performed,  with  each  eatiaation  cycle  aeparated  by  60 
aeconda  of  atate  propagation.  Initial  condltiona  were 
identical  to  Teat  Caae  1.  Reaulta  ahow  thia  type  of 
eatiaation  acheae  to  be  good  at  aaintaining  high  accuracy 
in  total  pointing  angle  (aee  Figure  22). 

A  final  Teat  Caae  waa  run  to  inveatigate  the  long  term 
capability  of  the  eatiaator.  Uaing  identical  condltiona  to 


taat  Cam*  4,  raault a  ahow  that  tha  praaanca  of  noiaa  doaa 
not  aarioualy  dagrada  tha  aatiaator’a  parforaanca.  (aaa 
Figura  23). 
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TOTAL  POINTING  ERROR  (rad*)  PROJECTION  OF  I  COMPONENT 


TIME  (  sec  ) 
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e 

19 _ '  12 


Figure  20.  Teet  Caee  #6  State  Propagation 


I 


Figure  21.  Teat  Caee  #6  State  Propagation 


Chapter  Five  Cone luaiona 


The  goal  of  this  reeearch  was  the  development  of  a  fine 
attitude  estiaation  routine  for  a  torque-free  axisymaetric 
spinning  satellite.  The  results  obtained  in  the  previous 
section  demonstrate  that  the  model  developed  in  this 
research,  when  run  with  short  propagation  cycles,  can 
adequately  predict  the  satellite  model  torque-free  attitude 
state. 

Figure  22  shows  that  with  estimation  updates  each  60 
seconds,  the  total  pointing  error  was  reduced  from  an 
initial  error  of  0.0107  radians  and  maintained  below  a 
f  maximum  total  pointing  error  of  0.0004  radians.  This 

result  equates  to  a  ground  distance  error  of  approximately 
nine  statute  alias  from  geosynchronous  altitude,  and  it 
should  prove  suitable  for  aost  pointing  requirements. 

Figure  23  shows  what  happens  when  the  estiaation  update 
time  is  extended  to  five  minutes.  At  the  end  of  the  the 
last  propagation  cycle  in  this  figure,  the  total  pointing 
error  has  grown  to  approximately  0.00544  radians,  or  about 
120  statute  ailes  ground  error  froa  geosynchronous 
altitude.  This  error  is  about  half  as  large  as  the  initial 
estimate  error,  and  it  aight  present  an  outer  bound  for 
estimate  update  cycle  tiaes.  The  aodel  was  not  tested  for 
large  Initial  estimate  errors,  and  further  investigation  in 


this  area  seems  warranted.  Also,  the  package  developed 
might  not  be  adequate  for  an  environment  where  continous 
torques  are  present,  but  small,  impulsive  torques,  such  as 
those  generated  by  fine  attitude  correction  Jets,  could  be 
handled  by  the  torque-free  model. 

The  attitude  prediction  capability  of  this  estimation 
package  was  shown  to  be  degraded  in  the  presence  of  a  low 
density  star  field.  This  problem  arises  from  the  method 
chosen  to  govern  the  estimation  sequence.  A  possible 
alternative  to  this  scheme  would  be  to  base  the  estimation 
sequence  on  a  predetermined  number  of  observations,  rather 
than  the  current  method,  which  searches  for  a  specified 
time  period  and  uses  all  data  recorded  during  that  period. 


Chapter  Six  Racoaaendatlona  for  Follow-On  Analyaia 

The  author  faala  thara  ara  aavaral  araaa  whara  tha  work 
praaantad  in  thia  paper  could  be  uaad  aa  a  baaia  for 
follow-on  reaearch.  Tha  firat  area  of  invaatigation  which 
praaanta  itaalf  would  be  tha  davalopaant  of  a  atar 
identification  aya tea.  Thia  addition  would  aatabliah  a 
coaplata  attitude  aatiaation  package. 

Another  area  for  follow-on  reaearch  aight  be  an 
invaatigation  of  alternative  aatiaation  routinea.  A a 
atated  earlier,  the  capability  of  thia  eatiaator  ia 
aerioualy  degraded  in  the  preaence  of  a  low  denaity  atar 
field.  Thia  drawback  could  be  overooae  by  the  choice  of  a 
different  baaia  for  the  leaat-aquarea  eatiaation  cycle. 
Another  alternative  aight  be  the  uae  of  a  Bayea  eatiaator 
inatead  of  the  leaat-aquarea  algoritha.  Still  another 
alternative  eight  be  the  uae  of  a  Kalaan  filter  to  replace 
or  augaent  the  batch  eatiaator. 

A  final  recoaaended  follow-on  reaearch  would  be  an 
invaatigation  of  the  uae  of  thia  aatiaation  package  a a  a 
controller  for  aoae  attitude  correction  package.  The 
reault  of  auch  an  effort  could  be  the  preaentation  of  an 
attitude  control  ayatea. 
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MAIN  PROGRAM  ESTMATOR. F 


PROGRAM  LISTING  eatmator.f 


program  estmatar 
This  is  tha  main  program  shall 
common/sidaraal /tgo 
common/dyncom/xgas0(3) , a(6) 

common/ ham/ t , xtruat 12, 4) , f ( 12, 4) , arrast ( 12 ) , n, h 
common/cov/q 

doubla  pracision  x0(6),x(6), t,h, tand, t0,xgas(8),xdal(6), tgo 
doubla  pracision  a, f , arrast , xtrua, xgasO, tsince, tnaw, tnow, q 
doubla  pracision  tsaad 
opan( 1 , f i la* ’ guassin *  ) 

raad  in  numbar  of  odas,no.  of  epochs,  no.  of  scans  between  update 
read( 1,104)  n, nepoch, nscan 

raad  in  start  time, timestep,  and  rel.  sidereal  time 
read( 1,105)  t,h,tgo 

specify  an  interval  to  be  scanned  and  initialize  time  parms 
tand* 12. 00d+00 

to*t 

taince*t 

tnow*t-tgo 

raad  initial  conditions 
raad( 1, 106)  ( x0( i ) , i*l , n ) 
write  out  initial  conditions 
writee  *, 1 ) 

writee*, 2)  t , ( i , x0( i ) , i*l, n ) 
raad  in  initial  guess 
raad( 1,100)  maxit 

100  format< lx, i2) 
do  30  i*l, 6 
read(l, 101)  xgas(i) 

30  continue 

101  formate  lx, e20. 13) 
raad  q  e sigma  squared) 
raade 1,101)  q 

raad  noise  flag  e if  flag*0  no  noise  will  be  added) 
raade 1 , 100 )  inoisa 

raad  noise  seed 

raade 1,101)  tsaad 
closae 1 ) 

writee *,107)  q,tgo 

if e inoisa. aq. 0)  writee * , 108 ) 

if e inoisa. na. 0)  writee * , 109 )  tsaad 

open  file  (errlist)  for  writing  errors 

opane  3, f 1 le* ‘errlist ’ ) 

do  300  iloop*l, nepoch 

iter*l 

do  46  ii*l , 6 
xgasOe ii )*xgase ii ) 

46  continue 

write  out  estimate  at  epoch 
writee  * , 40) 


40  foraatt 12x, ’ Initial  eatiaate  at  epoch  »*) 
write(*,41)  t, ( igea, xgea( igea) , iges=l,6) 

41  foraat(6x, ’tO  : ’ ,3x,e20. 13, /, 

+  6(8x, ’xgaa’, il, ’0  t  ’ , 3x, a20. 13, / ) ) 

writ#  out  invaatigation  tiaa  Info 
write(*,lll)  nepoch.nacan 

111  foraat(x, ///,6x, ’Thia  run  covara  ’,15,’  epoch  updatea’,/, 

+  6x, ’with  ’,15,’  acana  between  each  update’) 
call  caperr ( xO , xgea , tgo , t per r , r t go ) 
write(3,103)  rtgo.tperr 
write ( #, 4) 

call  truth  (tO, tend,h,xO, t,x) 
tnew»t 

now  order  truth  data  by  time 
write(*,5) 

call  order ( inoiae, taeed) 
get  tiae  realduala 
10  continue 
write(*,6) 
cell  reaid ( to, xgea) 
get  T 
write( *,7) 
call  dtdxo(tO.xgea) 
get  atate  correctiona 
write(*,8) 

call  covar(xdel, iter, if ail ) 
correct  the  gueaa 
do  9  i*l , 6 

xgea( i )«xgea( i )-xdel ( i ) 

9  continue 

write  out  updated  gueaa 

write<*,42)  tnow, ( J , xgea( j ) , xdel ( J ), J«l,6) 

42  foraat(2x, ’EST  Correctiona  and  Updated  ESTIMATE’,/, 

+  4x, ’tepoch  t  ’ , 3x, e20. 13, /, 

+  3(4x, ’xgea’ , 11, ’0  <  ’ , 3x, e20. 13, ’  froa  xdel  of  :’,3x, 

+  e20. 13, /  )  ) 

cell  caperr(xO, xgea, tgo, tperr , rtgo ) 
wrlte(3,103)  rtgo, tperr 
if  no  convergence  go  again 
itersiter'*'l 

if ( iter. le. aaxit. and. ifail . eq. 1 )  goto  10 
ifCifail.eq. 1)  write<  *, 102)  aaxit 
102  foraat(2x, ’PROGRAM  FAILED  AFTER ’, 12, ’ ITERATIONS  EXCEEDED’) 
At  thia  point  we  have  an  eatiaata  for  epoch  atate 
and  the  true  atate  at  tepoch  and  tend 
now  need  to  bring  the  eatiaate  forward  to  tend 
flrat  update  tiaea 
t»tnew-tend 
tgo*tgo-tend 
taince«taince+tend 

aet  eatiaated  atate  to  leaat  aquare  converged  value 
do  50  ip*l,6 


xgesO( ip)*xges( ip) 

50  continue 

cell  approx2( tsince, tO) 
write  out  finel  conditions 

write(*,3)  tnew, tsince, (i,x(i),a(i),i*l,6) 
cell  cmperr(x, a, tgo, tperr, rtgo ) 
write(3,103)  r tgo, tperr 

if  there  are  sore  then  one  scans  between  epoch  updates, 
perform  operations  to  bring  true  state  and  state  estimate 
forward  to  new  epoch 
if (nscan. eq. 1 )  goto  320 
inner  loop 

do  310  J loop*l , nscan-1 
initialize  halting 
nxt*0 

do  330  k»l,  6 
xtrue(k, 1 )*x(k) 

330  continue 

cell  heaing(nxt) 
if(nxt.ne.O)  goto  340 
write( * , 341 ) 

341  format(2x, ’HAMING  WOULD  NOT  INITIALIZE’) 
stop 

340  continue 
360  continue 

call  haaing(nxt) 

if ( t. le. tend)  goto  360 

update  times 

t*t-tend 

tgo*tgo-tend 

tsince* tsince-*- tend 

do  360  Ji»l,6 

x( Ji  )*xtrue< ji,nxt) 

350  continue 

bring  state  estimate  forward 
call  approx2( tsince, tO ) 
write  state  and  state  estimate 

write( * , 346)  tgo, tsince, ( xtrue( ii , nxt ),a(ii),ii-l,6) 

346  format(2x, ’State  propogation ’,/, 5x, ’sidereal  time  :  ’,e20.13, 

+  5x, ’time  since  epoch  :  ’, e20. 13, /, 20x, ’ STATES ’,/ , 

+  6(5x, e20. 13, 6x, e20. 13, / ) ) 
call  cmperr( x, a, tgo, tperr, rtgo  ) 
write(3,103)  rtgo, tperr 
310  continue 
320  continue 

now  reset  parameters  for  new  epoch 

do  43  i*l , 6 

xges< i )*a< i  ) 

x0( i )*xtrue( i, nxt ) 

43  continue 

to*t 

tsince»0. 00d+00 


tnow*t-tgo 
300  continue 

put  tflag  at  end  of  errlist  end  close  errlist 
rt «o *66400.0 
write! 3, 103)  rtgo,tperr 
103  formete2x,fl6.7,2x,fl6.7) 
close(3) 

now  cell  ploterr  which  puts  pointing  errors  into  errplot 
cell  plterr 

1  for«et(2x, ////,30x, 'MAIN  OUTPUT’,/) 

2  format e8x, ’ Initial  conditions  :’,/, 

+  8x, ’tO  : ’ ,3x, e20. 13, /, 

+  6(7x, ’x’.il, ’0  i ’,3x,e20.13,/)) 

3  format e8x, ’Final  conditions  :’,/, 

+  19x, ’true  state’ , 13x, ’state  estiaate’,/, 

+  8x, 't  t ’ ,3x,e20. 13, 'test: ' ,e20. 13, /, 

+  6(7x, ’x’.il,  ’  : ’ , 3x , e20 . 13 , Sx , e20 . 13 , / ) ) 

4  foraat( lx, ’ENTERING  TRUTH’) 


5  formate  lx, ’ENTERING  ORDER’) 

6  formate  lx, ’ENTERING  RESID’) 

7  foraet( lx, ’ENTERING  DTDXO ’ ) 

8  formate  lx, ’ENTERING  COVAR ’ ) 

104  foraate2x, il,2x, 13, 2x, 13) 

106  formate  lx, e20. 13, /, lx,e20. 13, /, lx,e20. 13) 

106  formate  lx, e20. 13, / , 4e lx, e20. 13, / ) , lx, e20. 13) 

107  formate 8x, ’covariance  :  ’, e20. 13, /, 8x, ’ initial  tgo 

106  formate 8x, ’NO  NOISE  ADDED’) 

109  formate8x, 'NOISE  rnd  SEED  :  ’,e20.13) 

stop 
end 


’ , e20. 13, /  ) 


include ’ truth ' 
inc 1 ude ’ torder ’ 
inc 1 ude ’ tstate ’ 
inc 1 ude ’ dtdxO ’ 
inc 1 ude ’ tresid ’ 
Include 'cover ian ’ 
inc 1 ude ’ coaperr ’ 
include 'ploterr ’ 


b* 
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PROGRAM  LISTING  coiparr 

subroutine  captrr(xtrut, xest, tgo, tperr , time ) 

this  routine  tekee  the  stete  vectors  (true  end  estimeted) 

end  computes  the  totel  pointing  error 

error  informetion  is  temporerily  held  in  file  (errlist) 

double  precision  xtrue(S) , xest(6) 
dimension  x(8),xg(8),err(3) 
do  1  i » 1 , 8 
x( i )»xtrue( i ) 
xg( i )*xest( i ) 

1  continue 
time*tgo 
do  2  i*l,3 

err(  i  )*x(  i+3  )-xg(  i+3 ) 

2  continue 

tperr «sqrt( err < 1 )**2+err(2)**2+err(3)**2) 

return 

end 
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PROGRAM  LISTING  covarian 

subroutine  covar(d«lx, JreJ , if si 1 ) 
c  This  routins  rssds  ths  residual  info  and  T  froa  dtfile  and 

c  after  computing  the  covariance,  computes  guess  corrections 

coamon/cov/q 

double  precision  q, tres, Tval(8),ttrue, tat, qinv, rejec, tqmag 
double  precision  pinv(6, 8 ) , p< 8, 8 ) , delx(8 ) , tqinvr(8 ) , tmat(8, 6 ) 
dimension  work<72) 
c  constants 

qinv«l . OOd+OO/q 
c  set  rejection  criterion 

c  if  this  is  less  than  the  third  iteration,  rejec  set  to 

c  3  seconds  (l£05*sigma) 

c  if  third  or  higher  iteration,  rejec  set  to  3*sigma 

c  if  all  residuals  are  rejected,  rejec  flag  is  widened 

rejec*3. 00d+00*daqrt( q ) 
if ( JreJ . 1 t. 3 )  rejec*3. 00d+00 
30  continue 

c  open  file  dtfile  for  reading 

open( 1, file* ’dtfile ’ ) 

c  initialize  tqinvr  vector  and  pinv  matrix 

do  2  i-1,8 
tqinvr( i )*0. 00d+00 
do  1  J  *1 , 6 
pinv( i, J ) *0. OOd+oo 

1  continue 

2  continue 
ifai 1 *0 

10  continue 

c  read  in  a  data  line 

read(l, 100)  isl it, ttrue, tst 

100  format(7x, 11 , 2(2x, e20. 13 ) ) 
if ( islit. eq. 0)  goto  20 
read(l, 101)  (Tval( i ) , i*l,6) 

101  format(3(2x, e20. 13 ) , / , 3(2x, e20. 13 ) ) 

c  compute  residuals  and  print  them 

tres*ttrue-tst 

irej*0 

c  reject  residual  if  limit  is  exceeded 

if (dabs(tres). gt. rejec)  then 
irej*l 

write(*,202)  tres 
tres*0. 00d+00 
endif 

if  (irej.eq.l)  goto  10 
write<*,102)  tres 

102  format( 2x, ’RESIDUAL  ’,e20,13) 

202  format(2x, ’RESIDUAL  ’,e20.13,’  This  residual  rejectee 
c  compute  tqinvr  and  pinv 

do  4  i*l,6 

tqinvr ( i ) * tqinvr ( i )+Tval ( i ) *trea*qinv 


07 


20 


o 

c 


c 

c 

c 

c 


do  3  J«l,6 

pinv( i, J )»Tval ( i )#Tval < J )*qinv+pinv( i, J ) 

continue 

continue 

now  get  enother  dete  line 
if ( islit. ne. 0)  goto  10 
continue 
closet  1 ) 

if  ell  residuels  heve  been  rejected  et  this  point 
need  to  widen  rejection  criterion  end  repeet  routine 
tq*eg*dsqrt ( tqinvr ( 1 ) *  *2+tqinvr ( 2 ) **2+tqinvr ( 3 )  *  *2 
+  +tqinvr <  4 ) *#2+tqinvr <  5  )  *  #2+tqinvr  (  8 )  *  *2 ) 

if ( tqseg . eq . 0 . 00d+00 )  re J ec * r e J  ec * 1 . 00d+0 1 

goto  30 

e  successful  residuels  pess 
inverting  pinv 

is  used 


299 

300 


c 

c 


c 

c 


c 

c 


103 


9 

104 


if ( tqeeg. eq. 0. 00d+00) 
et  this  point  we  heve 
coepute  coverience  by 
here  en  issld  routine 


edd  -liasld  to  f77  statement  to  link  it  to  the  program 

do  300  J*l«6 

do  299  i-1,6 

tmat( i, J )*pinv(i, J) 

continue 

continue 

call  leqt2f(pinv,6,6,6,p,7,work, ier) 
call  1 invlf ( tmat , 6 , 6 , p , 0 , work , ier ) 
call  ppinv(p.pinv) 

compute  state  corrections  at  epoch 
do  6  1*1,6 
delx( 1 )*0. 00d+00 
do  5  J  *1 , 6 

delx( i )*delx( i )+p( i, J )*tqinvr(J ) 

continue 

continue 

check  for  convergence 

Dr  Wiesels  fast  and  dirty  check 

do  7  i*l,6 

write  delx 

write(*,*)  delx(i) 

if (dabs(delx(l) ).gt.0. 10d+00*dsqrt< daba( p( i , i ) ) ) )  ifail=l 
continue 

if ( if ail . eq. 0)  then 
write  out  final  p  matrix 
write<*, 103) 
format( 2x, ’p  matrix’) 
do  9  i*l,6 

write( * , 104 )  < p< i , J ) , J *1 , 6  ) 
continue 

format(lx,6(lx,el2.5) ) 

end  if 

return 

end 

include 'debug ’ 


es 


a 


» 


$ 


a 


3 
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PROGRAM  LISTING  debug 

b  subroutine  ppinv(p, pinv) 

c  this  subroutine  multiplies  p  by  pinv  to  see  if  it  is  identity 

double  precision  p(6,8),pinv(8,6),prod(6,6),sum 
do  2  i*l,  8 
do  1  d*l,6 
sum»0. 00d+00 

i  do  5  k*l,6 

*um*su*+p( i , k)*pinv(k, j ) 

6  continue 

prod(i, j )*sua 

1  continue 

2  continue 

£  write<*,10) 

10  formet(2x, ’P  MATRIX’) 
writeC*,*)  ((p(i,d),d*l,S),i*l,6) 
write( *, 11 ) 

11  format (2x, ’PINVERSE  MATRIX’) 
write<*,*)  ( (pinvd,  j  ),J*1,S),  i*l,8) 

£  write(*,12) 

12  formet<2x, ’PRODUCT  P*PINV  (should  be  identity)’) 
write( *, * )  ( (prod( i, j ) , j=l, 8) , i»l. 6) 

return 

end 

* 


PROGRAM  LISTING  dtdxO 


subroutine  dtdxo(tO.xO) 

this  routine  finds  the  pertiels  of  the  estimate  time 
with  respect  to  the  initial  state 

double  precision  t0,x0(6),xpass0(6),xdel, ttrue, tstl, tst2 
double  precision  sdat(2),Tval(6),xst(6) 
del =0. 10d-02 

open  statfile  for  reading  and  dtfile  for  writing 
open( 1, file* ’statf ile ’ ) 
open(2, file* ’dtfile ’ ) 

10  continue 

get  ttrue, tat, sdattRA, dec ), islit  from  statfile 
read(l, 100)  islit, ttrue 

100  format(7x, il,4Sx,e20. 13) 

if  islit*0  we  have  reached  the  end  of  the  file 
if ( islit. eq. 0)  goto  999 
read(l, 101)  (sdat( i ) , i*l , 2 ) 

101  format(24x, 2(2x, e20. 13) ) 
readd,  102)  tstl 

102  formatdx,  /,5x,e20. 13,  /  ) 
do  2  idel =1,6 

store  originel  state  in  xpassO 
do  1  i=l , 6 
xpass0( i)«x0( i ) 

1  continue 

vary  the  initial  state  component 
xdel*x0( idel )*del 

if (x0< idel ) . eq . 0 . 00d+00 )  xdel=del*0. 10d+00 
xpass0( idel )=xpass0( idel )-xdel 

call  tstate( to, tstl, xpassO, sdat, islit, xst, tst2) 

Tval ( idel )  =  ( tst2-tstl )/xdel 

2  continue 

write  data  out  to  dtfile 

write( 2, 200 )  isl it , ttrue, tstl , ( Tval d ) , i*l , 6 ) 

200  format(2x, ’Slit  ’ , il , 2( 2x, e20. 13 ) , / , 

♦  3<2x,e20. 13), /,3(2x,e20. 13) ) 

get  another  siting  and  repeat  process 
if ( islit. ne. 0)  goto  10 
999  continue 

write  out  end  of  data  flag 
write(2,201)  isl it, ttrue, tstl 

201  format (7x, il , 2( 2x, e20. 13 ) , / , lOx, ’ END  OF  DATA’) 
closet  1 ) 

closed) 

return 


PROGRAM  LISTING  dynol 


*, 


subroutine  approx(t.tO) 
common  /dyncom/  z(6),a(6) 

c  thim  subroutine  computes  firmt  ordmr  dynamics  from; 

c  t  -  tima  (from  call) 
c  tO-  tima  at  stata  0  (from  call) 

c  z(i)  -  initial  stata  (from  common  aprox) 
c  outputs  are: 

c  a(i)  -  first  ordar  stata  approximation  at  tima  t 

c  NOTE  -  statas  1,2,3  ara  axact 
c 

doubla  pracision  ak,  am,  z,  a,  t,  to,  ta,  tl,  d 
c  dafina  constants: 

ak*-1.00d+00 
am*7.2921162d-06 

c  approximate  solutions  to  aquations  of  motion 

ta*t-tO 

a( 1 )*z(l )*dcos(z(3)*ak*ta)+z(2)*dsin(z(3)*ak*ta) 
a(2)«-z(l )*dsin(z(3)*ak*ta)+z(2)*dcos(z(3)*ak*ta) 
a(3)*z(3) 

d« 1 . 0d+00/ ( ( ak-1 . 00d+00 ) *z( 3 ) ) 
tl * ( ak-1 . 00d+00 ) *z( 3 ) * ta-z ( 6 ) 
a(8)»z(6)+z(3)*ta 

a(6)»z(6)+am*ta+d*(z( 1 )*(dcos(tl )-dcos(z(6) ) ) 

+  +z(2)*(dsin( tl )+dsin(z(6) ) ) ) 

a( 4 ) *z( 4 )+d* ( z( 1 )* ( dsin( tl )+dsin(z(6 ) ) ) 

+  -z(2)*(dcos( tl )-dcos(z(S) ) ) ) 

return 
and 


an 


V 


PROGRAM  LISTING  dyno2 


subroutine  approx2( t , tO ) 
common  /dyncom/  z(8),a(6) 

this  subroutine  computes  second  order  dynemics  from: 

t  -  time  (from  cell) 

to-  time  mt  stete  0  (from  cell) 

z(i)  -  initiel  stete  (from  common  mprox) 

outputs  are: 

m(i)  -  second  order  state  approximations  at  time  t 
NOTE  -  states  1,2,3  are  exact 

double  precision  z, a, ak, am, tl, tO, ta, t, d, a6, a4 

double  precision  ql , q2, q3, q4, q6, q6, q7 

double  precision  rl , r2, r3, r4, r5, r8, r7 

double  precision  si , s2, s3, s4, s5 

define  constants: 

ek*-l . 00d+00 

am*7. 2921152d-05 

approximate  solutions  to  equations  of  motion 
ta»t-tO 

a(  1  )*z(  1  )*dcos(z(3)*ek*ta)+z(2)*dsin(z(3)*ak*ta) 
a(2)»-z( 1 )*dsin(z(3 )*ak*ta)+z(2)*dcos(z(3 )*ak*ta ) 
a(3)*z(3) 

d*l .  0d+00/  ( ( ak-1 . 00d+00 ) *z( 3 ) ) 
tl»(ak-l. 00d+00)*z(3)*ta-z(8) 
a(6)*z(6)+z(3)*ta 

a8>z( 4 )+d* ( z( 1 ) *dsin(z(8) )+z( 2 ) *dcos( z( 6 )  )  ) 
a4«-z(4)*d*(z( 1 )*dcos(z(8) )-z( 2 ) *dsin( z( 6 ) ) )-d*d* 

♦  ( ( z( 1 ) **2-z( 2 )**2)*dsin(z(6 ) )*dcoa(z( 6 ) )/2. 00d+00+ 

♦  z( 1 )*z(2)*(dcos(z(6) ) )**2) 
sl«-(d*(z(l)**2+z(2)#*2)/2.00d+00)*ta 
s2*( a8*d*z( 1 ) )*( dcos( tl )-dcos( z( 6 )  ) ) 
s3»(a6*d*z(2))*(dain(tl  )+dain( z( 8 )  ) ) 
s4»( d*d*(z( 1 )**2-z( 2)**2 )/4. 00d+00)* 

+  (dsin(2.00d+00*tl)+dsin(2.00d+00*z(6) ) ) 

s5«(d*d*z(l )*z(2) )*( (dsin(tl ) )**2-<dain(z(S) ) )**2) 
a(6)sa(6)+sl+s2+a3+s4+s5 
end  of  a(8)  calculation 
start  a(5)  calculation 

a(5)*z(5)+aa*ta+d*(z( 1 )*(dcoa( tl )-dcos(z(6) ) ) 

♦  +z(2)*(dsin( tl )+dsin(z(6) ) ) ) 
rl»-(d*d*(z(l)**2+z(2)**2)/2.00d+00) 

+  * ( z( 1 )*dsin( tl )-z( 2 )*dcos( tl ) ) *ta 

r2*( ( a8*e6*d*z( 2 )/2. 00d+00 )+( a4*d*z( 1 ) >-(d**3*z( 1 )**2 
+  *z( 2 ) /2. 00d+00 ) ) * ( dsin( tl )+dain(z(6) ) ) 

r3* ( ( a6*a6*d*z( 1 )/2. 00d+00 )-( a4*d*z( 2 ) )-(d*  *3*3. 00d+00 
+  *z(l)*z(2)**2/2.00d+00) )*(dcos(tl )-dcos(z(6) ) ) 

r4«(2.00d+00*d*a6*z( l)*z(2))« 

♦  ( (dsin(tl) )**2-(dsin(z(6) ) )**2) 
r6»(a8*d*d*(z(l ) **2-z(2 ) **2 ) /2. 00d+00 ) * 

+  (dsin(2.00d+00*tl)+dain(2.00d+00*z<6) ) ) 


r8*d*d*d*(z(l)**2*z(2)-(z(2)**3/3. 00d+00) )* 

(  (dain(tl)  )**3-Md*in<z(6) )  )**3> 
r7*d*d*d*(z<l)*z(2)**2-(z(l)**3/3.00d+00) )* 

( <dcoa< tl ) )**3-<dcoa(z(6) ) )**3) 
a(5 )*a(5 )+rl+r2+r3+r4+r6+r6+r7 
•nd  atata  5  calculation 
bag in  atata  4  calculation 
a(4)*z(4)+d*(z(l)*(dain(tl)+dain(z(6) ) ) 

-z( 2 ) * ( dcoa( tl )-dcoa(z(6) ) ) ) 
ql * ( d  *d* ( z ( 1 ) *  * 2+z ( 2 ) *  *2 ) / 2 . 00d+00 ) * 

(z( 1 ) *dcoa( tl )+z(2)*dain( tl ) )*ta 
q2»-(d*a4*z(2)+(d*d*d*(z(l)**2+z(2)**2)*z(l >/2. 00d+00) )* 
(dain( tl )+dain(z(6) ) ) 

q3*-(d*a4«z< 1 )+<d*d»d*(z( 1 )**2-z( 2)**2)*z(2)/2. 00d+00 )  )* 
(dcoa( tl )-dcoa(z<6) ) ) 
q4* ( d*d*aS* (z( 1 )#*2-z<  2 )**2 )/2. 00d+00  )* 

( (dain(tl ) )**2-<dain(z(6) ) )**2) 
q5*-( d*d*a6*z<  1  )*z(2)/2.  00d+00)* 

( dsin( 2. 00d+00*tl >+dain< 2. 00d+00*z( 6 ) ) ) 
q6«d*d*d*(z( 1 )*»3/fl. 00d+00-z( 1 )**2*z(2)/2. 00d+00 )* 

( (d«in( tl ) )**3+(dain(z(6) ) )#*3) 
q7»-(d*d*d*(z(  1 )  **2+z(  2 )  **Z )  *z(  2 ) /6.  00d+00 )  * 

( (dcoa(tl ) )**3-(dcoa(z(6) ) )**3) 
a(  4  )  »a(  4  )-*-ql+q2>q3+q4-*-q5+q6+q7 
raturn 
and 


PROGRAM  LISTING  gutsaln 


8  5  30  no.  odes, epoch  update*, scene  between  updat 

0. 0000000000000d+00  start  tise 

0. 1000000000000d+00  nuaerical  integration  tiaestep 

0. 0000000000000d+00  initial  sidereal  tiae 

0. 1000000000000d-01  initial  true  oaegalO 

0. 5000000000000d-01  initiel  true  oaeg*20 

0. 6223598778000d+00  initial  true  oaega30 

0. 5000000000000d-01  initial  true  psilO 

0. 5000000000000d-01  initial  true  psi20 

0.8000000000000d+00  initial  true  psi30 

10  aax  nuaber  of  residual  passes 

0.0110000000000d+00  initial  est.  oaegalO 

0. 0650000000000d+00  initial  est.  oaega20 

5.7236987780000d-01  initial  est.  oaega30 

0. 0580000000000d+00  initial  est.  psilO 

0 . 0660000000000d+00  initial  est.  psi20 

0.8080000000000d+00  initial  est.  psi30 

1 . 0000000000000d-09  covariance 

0  noise  flag  (0«no  noise) 

1 . 0000000000000d+00  noise  rnd  seed 


PROGRAM  LISTING  haming 


THIS  NUMERICAL  SUBROUTINE  PROVIDED  BY  DR  WILLIAM  WIESEL 


aubroutina  haming(nxt) 

haming  ia  a  fourth  ordar  predictor-corractor  algorithm 
for  tha  intagratlon  of  ayatama  of  ordinary  diffarantial  equations 
tha  common  /ham/  contalna  moat  of  tha  variables: 
x  ia  the  indapandant  variabla,  the  ’time* 
y  contalna  4  copies  of  tha  atate  vector,  with 
n  odea  being  integrated 
f  contalna  the  calculated  aquations  of  motion 
arrest  is  a  truncation  error  estimate 
n  ia  tha  number  of  ode  a 
h  is  the  integration  tlmestep 
nxt  assumes  tha  values  1, 2, 3,4, 1, 2, 3, 4. . . . ,  and  points  to 
tha  currant  value  of  tha  state  vector 

tha  user  must  supply  a  aubroutina  ’rhs(nxt)’  which 
calculates  tha  aquations  of  motion  f(i,nxt)  from  tha 
state  vector  y(i,nxt) 

to  initialize  haming,  tha  initial  conditions  must  be  stored 
in  x(  1 , 1 ) ,  1*1, n  »  x,n,and  h  must  be  initialized,  and 
than  haming  is  called  with  nxt*0.  If  haming  returns  with 
nxt*l,  initialization  ia  successful.  If  nxt*0  still,  haming 
did  not  initialize  (h  is  usually  too  big) 

common  /ham/  x, y ( 12, 4 ) , f ( 12, 4 ) , errest( 12 ) , n, h 
double  precision  x, y, f, arrest, h, xo, tol , hh 
to 1  *  1.0d-12 

branch  on  nxt:  startup  or  propagating? 
if (nxt)  190,10,200 

haming  initialization:  4  point  picard  iteration 

10  xo  *  x 

hh  *  h/2. d+00 
call  rhs(l) 
do  40  1  *  2,4 
x  *  x  ♦  hh 
do  20  1  *  l,n 

20  y(i,l)  *  y( i, 1-1 )  +  hh*f(i,l-l) 
call  rha< 1 ) 
x  *  x  ♦  hh 
do  30  i  *  1 , n 

30  y(i,l)  *  y ( i , 1-1 )  +  h#f ( i, 1 ) 

40  call  rhs( 1 ) 


do  120  1  *  l,n 

hh  -  y( 1, 1 )  +  h*(  9. d+00*f ( 1 , 1 )  +  19. d+00#f ( 1, 2)  -  5. d+00*£( 1, 3) 
1  +  f(i,4)  )  /  24. d+00 

lf(  daba(  hh  -  y(i,2) )  .It.  fcol  )  go  to  70 
law  *  0 

70  y(i,2)  »  hh 

hh  *  y( 1, 1 )  +  h*(  f(i,l)  +  4.d+00*f(i,2)  +  f ( 1 , 3 ) )/3. d+00 
If (  daba(  hh-y( 1, 3) )  .It.  tol  )  go  to  90 
law  *  0 

90  y(i,3)  =  hh 

hh  «  y ( i , 1 )  +  h*<  3.d+00*f( 1,1)  +  9. d+00*£ ( 1 , 2 )  ♦  9. d+00*£( 1 , 3 ) 

1  +  3. d+00*f ( 1,4)  )  /  8. d+00 

If (  daba(hh-y(i,4) )  .It.  tol  )  go  to  110 
law  «  0 

110  y(i,4)  *  hh 
120  contlnua 
X  3  xo 

do  130  1  *2,4 
x  *  x  +  h 
130  call  rha( 1 ) 

if (law)  140,140,150 
140  Jaw  *  Jaw  +  1 

if (Jaw)  50,280,280 
150  X  *  XO 
law  *  1 

Jaw  *  1 

do  160  1  *  i,n 
160  arraat(i)  *  0.0 
nxt  ■  1 

go  to  280 

190  Jaw  -  2 

nxt  »  laba(nxt) 

haming  propagation  saction 

200  x  *  x  +  h 

npl  *  «od(nxt,4)  +  1 
go  to  (210,230), law 
210  go  to  (270, 270, 270, 220), nxt 
220  law  *  2 

parauta  indlcaa 
230  na2  *  «od(npl,4)  +  1 

nal  *  «od(na2,4)  +  1 

npo  *  aod(nal,4)  +  1 

pradictor 
do  240  1  *  i,n 

f(l,na2)  ■  y(i,npl)  +  4.d+00*h*(  2. d+00*f ( 1 , npo )  -  f(i,nal) 

1  +  2. d+00*f ( 1 , n»2 )  )  /  3. d+00 

240  y( 1 , npl )  *  f(i,n*2)  -  0. 925619835d+00*arraat ( 1 ) 
call  rha(npl) 


corrector 
do  250  i  *  l,n 

y(i,npl)  *  <  9. d+00*y< i, npo )  -  y(i,n*2)  +  3.d+00*h*<  f(i,npl) 
1  +  2. d+00*£( 1 , npo )  -  £<i,nal>  )  )  /  8.d+00 

errest(i)  *  f(i,n«2)  -  y(i,npl) 

250  y(i,npl)  *  y(i.npl)  +  0. 0743801663d+00  *  errest(i) 
go  to  (260*270 ) , jew 
260  cell  rhe(npl) 

270  nxt  ■  npl 
280  return 
end 


Include  ’theslrhs’ 


PROGRAM  LISTING  matxmat 


subroutine  matmat( id, rinl, rin2, rout ) 

This  routine  multiplies  two  square  matrices  together 
double  precision  rinl( 3, 3 ) , rin2( 3, 3 ) , rout( 3, 3 ) , sum 
do  15  i=l , id 
do  10  j  =1, id 
sum*0. 00d+00 
do  5  k*l,id 

sum*sum+r ini ( i , k ) *rin2( k, j ) 

6  continue 

rout( i, j )saum 
10  continue 
15  continue 
return 
end 
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PROGRAM  LISTING  noise 


subroutine  nois( tseed, t ) 

c  this  routine  adds  pseudo  randos  noise  to  the 

c  true  slit  crossing  time 

c  slit  crossing  tiae  and  noise  to  be  added  are 

c  written  to  an  output  file 

c  and  noisy  tiaes  are  passed  along 

coaaon/cov/q 

double  precision  tseed, t, q, tnoise 
c  pseudo  randoa  noraal  (0,1)  inforaation  is 

c  obtained  through  IMSL  routine  ggnqf 

rnoise*  ggnqf (tseed) 
tnoise*rnoise 
tnoise-tnoise*dsqrt ( q ) 
write(*,101)  t, tnoise 

101  foraat(5x, ’TIME  e20. 13, 2x, ’NOISE  :’,e20.13) 
t*t+tnoise 
return 
end 


PROGRAM  LISTING  observ2 


subroutine  obaer( x, t , n, sdat , obs ) 

This  routine  computes  the  values  of  the  observation 
relations  from  a  given  state, time, and  star  data 
common/sidereal /tgo 
double  precision  x(6),t,tgo 

double  precision  p(3),romeg(3,3),rl(3,3),r2(3,3),r3(3,3) 
double  precision  rthet(3,3),obs(2,3),sdat(2) 
double  precision  tmat(3,3), tvec(3),romega, thetaz, thetai 
double  precision  tmatl(3,3), tvecl(3) 

double  precision  hrl(3, 3 ) , hr2( 3, 3 ) , hr3(3, 3 ) , romdot(3, 3 ) , pro< 3 ) 
constants 

tgo  specifies  the  starting  orbital  position  of  the  satellite 
in  order  for  the  star  data  table  to  be  valid  in  this  case 
tgo  should  be  close  to  0 

thetai  and  thetaz  are  star  sensor  slit  geometry  parameters 

for  this  simulation 

thetai *15  degrees 

thetaz*0.5  degrees 

thetai*2. 617993878d-01 

thetaz*8. 72664626d-03 

romega  is  the  satellite  nominal  orbital  rate 
for  this  simulation 
romega3 1  rev  per  24  hrs 
romega37. 2921152d-05 

nf lag  (n)  from  call  tells  what  is  desired  from  this  routine 
n=l  observation  state  desired 
n=2  partial  obs  state  wrt  state  is  desired 

n=3  partial  obs  state  wrt  t  is  desired 

component  setup 
star  vector 

p( 1 )=dcos(sdat( 1 ) ) *dcos( sdat ( 2 ) ) 
p(2)=dsin(sdat( 1 ) ) *dcos( sdat ( 2 )  ) 
p( 3)=dsin(sdat(2) ) 
romeg  matrix 

romegd,  1 )  =  -dsin( romega* ( t-tgo ) ) 

romeg ( 1 , 2 )30. 00d+00 

romeg ( 1 , 3 )*-dcoa( romega* ( t-tgo ) ) 

romeg(2, 1 ) 3 -romeg ( 1, 3) 

r omeg ( 2 , 2 ) 3 0 . 00d+00 

romeg(2, 3)=romeg( 1,1) 

romegO,  1  )*0.  00d+00 

romeg(3, 2 )»-l . 00d+00 

romeg ( 3 , 3 ) 3  0 . 00d+00 

rl  matrix 

rl< 1 , 1 )»dcos(x( 5) ) 

r 1 ( 1 , 2 ) =0 . 00d+00 

rl( 1 , 3 ) 3dsin( x( 5) ) 

rl(2, 1 )30. 00d+00 

rl(2, 2)31. 00d+00 


r 1 ( 2 , 3 ) -0 . 00d+00 

rl(3, l)»-rl(l,3) 

rl(3,2)«0. 00d+00 

rl(3,3)*rl( 1, 1 ) 

r2  matrix 

r2( 1, 1 )»1. 00d+00 

r2( 1 , 2 )-0. 00d+00 

r2( 1 , 3 )-0. 00d+00 

r2(2, 1  )-0. 00d+00 

r2(2, 2)-dcoa(x(4) ) 

r2(2,3)«-dain(x(4) ) 

r2(3, 1 )*0. 00d+00 

r2(3,2)— r2(2,3) 

r2<3, 3)«r2(2, 2) 

r3  matrix 

r3( 1, 1 )*dco«(x(6) ) 

r3(l,2)*-dain(x(6) ) 

r 3 ( 1 , 3 ) *0 . 00d+00 

r3(2, 1 )*-r3( 1,2) 

r3(2,2)*r3( 1,1) 

r3(2, 3)-0. 00d+00 

r3(3, 1 )»0. 00d+00 

r3(3, 2)-0. 00d+00 

r3 ( 3 , 3 ) * 1 . 00d+00 

rthat  matrix 

rthat (1,1) -dcoa ( thataz ) 

rthat ( 1,2) >dcoa( thatai ) *dmin( thataz ) 

rthat ( 1 , 3 ) *-dain( thatai ) *dain( thataz ) 

rthat (2, 1 )--dain( thataz) 

rthat ( 2, 2 ) -dcoa( thatai ) *dcoa( thataz ) 

rthat ( 2, 3 ) *-dain<  thatai ) *dcoa< thataz ) 

rthat ( 3 , 1 ) -0 . 00d+00 

rthat (3, 2)-dain( thatai ) 

rthat (3, 3 )*dcom< thatai ) 

computa  obaarvation  matrix 

coaputa  obaarvation  ralation  1 

call  aatmat<3, r2, r3, tmat ) 

call  aataat(3, rl, tmat, tmatl ) 

call  aataat(3, romeg, tmatl , tmat ) 

do  3  ik-1,3 

tvec( ik)-p(ik) 

3  continue 

call  vecmat(3, tvac, tmat, tvacl ) 

do  2  J  -1 , 3 

oba( 1 , J ) -tvacl ( J  ) 

2  continue 

computa  obaarvation  ralation  2 
call  vecmat(3, tvacl, rthat, tvac ) 
do  6  .j»l,3 
oba(2, j )-tvac(J ) 

6  continue 

if  n-1  G  haa  bean  found  ao  return 


if(n.eq.l)  return 

if  n* 2  compute  H  matrix 

if(n.mq.2)  then 

Sine*  thm  lmft  half,  of  thim  matrix  is  zero 

only  tha  right  half  will  be  computed 
compute  h  matrix 
hrl  matrix 
hrl< 1, 1 )*rl(3, 1 ) 
hrl(l,2)*0. 00d+00 
hrl<l,3)-rl<l, 1) 
hrl (2, 1 )*0.00d+00 
hrl(2, 2)*0. 00d+00 
hrl(2,3)*0. 00d+00 
hrl (3, 1 )*-hrl< 1, 3) 
hrl(3,2)*0. 00d+00 
hrl(3,3)*hrl(l, 1) 
hr2  matrix 
hr2( 1, 1  )*0. 00d+00 
hr2( 1 , 2 ) *0. 00d+00 
hr2( 1 , 3 )*0. 00d+00 
hr2(2, 1  )*0. 00d+00 
hr2(2, 2 ) *r2(2, 3 ) 
hr2(2, 3 )*-r2(2, 2) 
hr2(3, 1  )*0. 00d+00 
hr2(3, 2)*r2(2, 2) 
hr2(3, 3)*r2(2,3) 
hr3  matrix 
hr3(l,l)*r3(l,2) 
hr3( 1, 2)*-r3< 1, 1 ) 
hr3( 1 , 3)*0. 00d+00 
hr3(2, 1 )*r3( 1, 1 ) 
hr3( 2, 2)*hr3( 1, 1 ) 
hr3 ( 2 , 3 ) *0 . 00d+00 
hr3(3, 1  )*0. 00d+00 
hr 3 ( 3, 2 ) *0. 00d+00 
hr 3 O , 3 ) *0 . 00d+00 
compute  h  matrix 
call  vecmatO,  p,  romeg,  tvec ) 
do  13  1*1,3 
pro( i )*tvec< i ) 

13  continue 

call  matmatO,  rl,  hr2,  tmat) 

call  matmatO,  tmat,  r3,  tmatl ) 

call  vecmat( 3, pro, tmatl, tvec) 

oba( 1, 1 )*tvec(2) 

call  vecmat(3, tvec, rthet, tvecl ) 

obs(2, 1 )*tvecl(2) 

call  matmatO,  hrl,  r2,  tmat) 

call  matmatO,  tmat,  r3,  tmatl ) 

call  vecmatO,  pro,  tmatl,  tvec) 

obs (1,2)*  tvec ( 2 ) 

call  vecmatO,  tvec,  rthet,  tvecl ) 


obs( 2 , 2 ) ■ tvacl ( 2 ) 

call  aataat( 3, rl , r2, taat ) 

call  aataat(3, tmat, hr3, taatl ) 

call  vacaat(3, pro, taatl, tvac) 

oba( l,3)“tvac(2) 

call  vacaat( 3, tvac, r that, tvac 1 ) 

oba(2, 3)*tvacl(2) 

now  Haatrix  has  baan  coaputad 

andif 

if  n*2  wa  ara  dona  so  raturn 

if(n.aq.2)  raturn 

if  n*3  than  coaputa  tfdot 

gdot  will  ba  passad  back  as  a  aatrix  but  only  row  2 
is  aeaninfful 
if(n.aq.3)  than 
roadot  aatrix 

roadot ( 1,1) sroaaga*roaef ( 1,3) 

roadot ( 1 , 2 ) *0 . 00d+00 

roadot (1, 3 )*roaafa*roaaf (2, 1 ) 

roadot ( 2 , 1 ) * r oaaga *  roaaf (1,1) 

roadot ( 2 , 2 ) *0 . 00d+00 

roadot (2,3) “roadot (1,1) 

roadot ( 3, 1 ) *0. 00d+00 

roadot ( 3 , 2 ) *0 . 00d+00 

roadot ( 3 , 3 ) *0 . 00d+00 

now  fat  tha  fdot  aatrix 

call  aataat(3,r2,r3, taat) 

call  aataat(3, rl, taat, taatl ) 

call  aataat( 3, roadot, taatl, taat) 

call  vacaat(3, p, taat, tvac) 

do  21  i*l,  3 

obs(l, i)stvac(i) 

21  continua 

call  vacaat(3, tvac, rthat, tvacl ) 

do  22  1*1,3 

oba( 2 , i ) “tvacl ( i ) 

22  continua 
andif 
raturn 
and 


Inc 1 uda ’ aatxaat ’ 
inc 1 uda ' vacxaat ' 


I 
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PROGRAM  LISTING  ploterr 


aubroutine  plterr 

c  thia  routina  reeda  fila  (arrliat)  and 

c  writaa  tha  inforaation  in  plottabla  format  in  file  (errplot) 

open(l,file« ’ arrliat ’ ) 
opan(2, f ila= ’ arrplot ’ ) 

c  acan  arrliat  for  ain/aax  error,  ain/max  t 
tain*0. 0 
taax»0.0 
errain*0. 0 
arraax*0. 0 
10  continue 

raad( 1,100)  t,tparr 
if (t. ge. 86399. 0 )  goto  20 
if ( t. It. tain)  tain*t 
if ( t. gt. taax)  taax=t 
if ( tparr. 1 t. errain)  errain*tperr 
if ( tparr. gt. arraax)  erraax=tperr 
if(t. It. 86399.0)  goto  10 
20  continue 
cloae( 1 ) 

c  aat  acala  factora 
tacale*taax-tain 
aacalaaarraax-arrain 
writa( *,101)  tacala, aacala 

c  read  fila  again  and  write  to  arrplot 

open( 1, f ila* 'arrliat ' ) 

30  continue 

raad< 1,100)  t, tparr 
if(t.ge. 86399.0)  goto  40 
tp lot* ( taax-t )/ tacala 
ap 1 ot * tperr /aacala 
write(2,102)  tplot,eplot 
if (t. It. 86399.0)  goto  30 
40  continue 
cloaa( 1 ) 
cloae(2) 

100  format (2x, fl5.7,2x,fl5.7) 

101  format(2x, /,5x, ’ERRPLOT  PARAMETERS',/, 

+  5x, ’Tima  acala  factor  :  ’,fl5.7,/, 

+  5x, ’Error  acala  factor:  ’,fl5.7,/) 

102  format (2x, 2f8. 5, "\”+\nn ) 
return 

and 


PROGRAM  LISTING  stardat 

1 

0. 1445190798496@+01 

-0.54066725450770-02 

2.23 

1852 

36 

2 

0. 14004958252700+01 

-0.69279875039600-02 

4.73 

1766 

26 

3 

0. 17141096762790+01 

-0.21079698867400-01 

5.  10 

2395 

107 

4 

0. 1 4635749333000+01 

-0.21137876499140-01 

1.70 

1903 

46 

5 

0. 148342078 131 la+01 

-0. 34033920418330-01 

1.77 

1948 

50 

6 

0. 14123931528960+01 

-0.42072131252180-01 

3.36 

1788 

28 

7 

0. 14746577740230+01 

-0.45519156525320-01 

3.  81 

1931 

48 

8 

0. 1689449627358a+01 

-0. 82932228301430-01 

5.06 

2344 

102 

9 

0.  13399183567380+01 

-0. 89113602736380-01 

2.79 

1666 

10 

10 

0. 1460287896624e+01 

-0. 10330894732110+00 

2.76 

1899 

44 

11 

0. 163231 1910943©+01 

-0. 10941759970390+00 

3.  98 

2227 

87 

12 

0. 1382533478382a+01 

-0. 11973443483930+00 

3.  60 

1736 

21 

13 

0. 16535540223760+01 

-0. 13641202547160+00 

5.27 

2273 

93 

14 

0. 14102260357960+01 

-0. 13651383634460+00 

4.  14 

1784 

27 

15 

0. 13691744373420+01 

-0. 14344667198540+00 

0.  12 

1713 

19 

16 

0. 13456633979290+01 

-0. 15312355306160+00 

4.27 

1679 

11 

17 

0. 15141603927690+01 

-0. 16885090887890+00 

2.06 

2004 

57 

18 

0. 16731089821790+01 

-0. 20108132240320+00 

5.22 

2305 

99 

19 

0. 13594951323090+01 

-0.20746631858430+00 

4.45 

1696 

13 

20 

0. 18041032148460+01 

-0 . 20976433643310+00 

4.07 

2574 

131 

21 

0. 13636039281180+01 

-0. 22618012667760+00 

4.36 

1705 

17 

2 

2 

0. 13912892132980+01 

-0.23024286532580+00 

4.29 

1756 

24 

2 

3 

0. 15520267664440+01 

-0.24731315603990+00 

3.71 

2085 

67 

2 

4 

0. 15108079061880+01 

-0.25877899859960+00 

3.  55 

1998 

54 

2 

5 

0. 13623749254640+01 

-0.28315058235220+00 

3.31 

1702 

14 

2 

6 

0. 1589529527668O+01 

-0.28767389399750+00 

4.93 

2148 

78 

2 

7 

0. 17648114899720+01 

-0.29136817424810+00 

-1.46 

2491 

119 

2 

8 

0. 1812720778029e+01 

-0.29728774929620+00 

4.38 

2596 

133 

2 

9 

0. 1448819628889O+01 

-0. 31124068703930+00 

2.  58 

1865 

39 

3 

0 

0. 16668621579240+01 

-0.31324296754260+00 

1.98 

2294 

96 

3 

1 

0. 17331483084830+01 

-0. 3180571673966e+00 

4.  43 

2443 

112 

3 

2 

0. 17279050485210+01 

-0. 33684013322200+00 

3.96 

2429 

110 

3 

3 

0. 1429337391052e+01 

-0.36252427823380+00 

2.  84 

1829 

32 

3 

4 

0. 15300210722890+01 

-0.36441989972720+00 

3.  81 

2035 

62 

3 

5 

0. 13299554346590+01 

-0. 39079861211980+00 

3.  19 

1654 

7 

3 

6 

0. 16070628144760+01 

-0.39137069226360+00 

5.50 

2180 

80 

3 

7 

0. 15001832143370+01 

-0.39187489849200+00 

3.60 

1983 

51 

3 

8 

0.  17209237314840+01 

-0.40058700034270+00 

4.54 

2414 

108 

3 

9 

0. 17069756419320+01 

-0. 40852340030350+00 

4.  34 

2387 

105 

4 

0 

0. 18041832089930+01 

-0.42173942125230+00 

3.86 

2580 

132 

4 

1 

0. 18239417906240+01 

-0. 50527766665510+00 

1.50 

2618 

136 

4 

2 

0. 16568337868180+01 

-0. 52457325116580+00 

3.02 

2282 

94 

4 

3 

0. 17857408966030+01 

-0. 56705747404700+00 

3.  96 

2538 

125 

4 

4 

0. 16911949564710+01 

-0. 56844888931200+00 

4.48 

2361 

103 

4 

5 

0.  16648113959980+01 

-0.58342963206020+00 

3.  85 

2296 

95 

4 

6 

0. 14795446959160+01 

-0. 59484214611500+00 

2.  64 

1956 

49 

4 

7 

0. 13828461831510+01 

-0. 60929444195080+00 

4.83 

1743 

22 

4 

8 

0. 16406094970890+01 

-0.61322143276830+00 

4.  37 

2256 

90 

4 

9 

0. 1557648180022©+01 

-0 . 61582973037300+00 

4.36 

2106 

68 

5 

0 

0. 14427836986100+01 

-0. 61926705937250+00 

3.  87 

1862 

35 

0.15289666025330+01  -o. 
0. 1775108932588@+01  -0. 
0. 1565000379366a+01  -0. 
0.17334901021830+01  -0. 
0.14387258081540+01  -0. 
0. 1328908236997a+01  -0. 
0.17870062602130+01  -0. 
0.13918419010610+01  -0. 
0.15137095160420+01  -0. 
0 . 16296047458300+01  -0 . 
0.16738071139070+01  -0. 
0. 17219127514220+01  -0. 
0. 1786860816248a+01  -0. 
0.  16144150138400+01  -0. 
0. 1525177783726a+01  -0. 
0. 1331875296809a+01  -0. 
0. 1780381281371a+01  -0. 
0. 1455 124630781 •♦01  -0. 
0. 15042338326990+01  -0. 
0.13690653642640+01  -0. 
0. 1906081348557a+01  -0. 
0. 1609360831213a+01  -0. 
0. 18711674914180+01  -0. 
0. 1796074700175a+01  -0. 
0.20159352805760+01  -0. 
0. 99252511256300+00  -0. 
0. 1617505701067a+01  -0. 
0. 12899074004700+01  -0. 
0. 14507758521 17a+01  -0. 
0 . 2 177068259431 •-*■01  -0 . 
0.2186584848196a+01  -0. 
0 . 22771 38347691a+01  -0 . 
0 . 6631305772321a+00  -0 . 
0. 28171577626830+01  -0. 
0 . 23524929380720+01  -0 . 
0. 5472610770051a+01  -0. 
0 . 49127600513050+01  -0 . 
0. 6135370464754e-*-01  -0. 
0 . 337388 1433046«-*-01  -0 . 
0. 59539144002130+01  -0. 
0 . 68396607840120+01  -0 . 
0.3868409783773a+01  -0. 
0.4324242299709a+01  -0. 
0 . 4267308205123a+01  -0 . 
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0. 15846644223360+01 

0. 40605569866630+00 

4.  16 

2134 

77 

271 

0. 16668912467440+01 

0.39309178083170+00 

2.88 

2286 

97 

272 

0. 16316283236450+01 

0. 39291239976970+00 

3.  28 

2216 

86 

273 

0. 13184580781820+01 

0.37644327902030+00 

4.64 

1602 

4 

274 

0. 14692108922040+01 

0.36886079304670+00 

3.  00 

1910 

47 

275 

0. 15422747380810+01 

0.35385096147760+00 

4.  41 

2047 

64 

276 

0. 16931584519840+01 

0.35295405616750+00 

4.  15 

2343 

104 

277 

0. 15838790242970+01 

0. 35150446326070+00 

4.63 

2135 

76 

278 

0. 14455834976190+01 

0.3243452008414O+00 

4.38 

1845 

37 

279 

0. 14235923489990+01 

0.31328175263710+00 

5.  42 

1808 

31 

280 

0. 15120296366090+01 

0 . 3093450655460e+00 

5.49 

1990 

55 

281 

0. 17314393402670+01 

0.28646670793140+00 

1.93 

2421 

111 

282 

0. 13250612406030+01 

0. 26848981663350+00 

4.68 

1638 

5 

283 

0. 15999724143880+01 

0.25780452310050+00 

4.42 

2159 

79 

284 

0. 16190474086190+01 

0.24806461624670+00 

4.  48 

2199 

84 

285 

0. 17589210037810+01 

0.23115916318320+00 

4.  49 

2478 

116 

286 

0. 17646151405150+01 

0.22537533696680+00 

3.36 

2484 

118 

287 

0. 16387478126640+01 

0.2142827989416O+00 

5.04 

2241 

89 

288 

0. 14585862004930+01 

0. 17321908014620+00 

3.  39 

1879 

43 

289 

0. 1746874667614O+01 

0. 1729766733067O+00 

4.65 

2456 

114 

290 

0. 13469688304660+01 

0. 17122649591660+00 

5.  43 

1672 

12 

291 

0. 15774722113340+01 

0. 16839518401860+00 

4.  12 

2124 

75 

292 

0. 14572117536960+01 

0. 16545721311070+00 

4.  41 

1876 

42 

293 

0. 14663092824370+01 

0. 16201988411120+00 

4.09 

1907 

46 

294 

0. 15460635570270+01 

0. 12924163112710+00 

0.50 

2061 

66 

295 

0. 17107062831880+01 

0. 12819928171250+00 

4.  50 

2385 

106 

296 

0. 14150184190060+01 

0.  11059084881230+00 

1.64 

1790 

29 

297 

0. 1439700283487e+01 

0. 10362407621390+00 

4.20 

1839 

34 

298 

0. 16709200483800+01 

0. 8031423442309e-01 

4.  33 

2298 

98 

299 

0. 13634512117810+01 

0.49630376541660-01 

4.  46 

1698 

16 

300 

0. 17760979524910+01 

0.42411500829000-01 

4.47 

2506 

123 

301 

0. 15621933082310+01 

0.96380959817160-02 

5.22 

2103 

69 

999 

0.  0+00 

0.  @+00 

0. 

0 

0 
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PROGRAM  LISTING  theslrhs 
subroutine  rhs(nxt) 

common  /hem/  t , x( 12, 4 ) , f { 12, 4 ) , err ( 12 ) , n, h 
double  precision  t, x, f , err, h, ak, aomega 
c  it  is  the  Job  of  rhs  to  celculete  the  current  equations 

c  of  motion  f(i,nxt),  1*1,  n  from  the  current  value  of 

c  the  state  vector  x(i,nxt),  i  =  1 , n,  at  the  current  time 

c  t.  The  other  three  copies  of  the  state  vector  are  ignored. 

ak*-l. 00d+00 
aomega*7. 2921152d-05 
f ( 1, nxt )=ak*x(2, nxt )*x(3, nxt ) 
f ( 2, nxt ) =-ak*x( 1 , nxt ) *x( 3, nxt ) 
f (3,nxt)=0.00d+00 

f ( 4, nxt ) sx( 1 , nxt ) *dcos( x(8,nxt))-x(2, nxt )*dsin(x(6, nxt ) ) 
f ( 5, nxt )*x( 1 , nxt ) *dsin(x(6, nxt ) ) /dcos( x( 4, nxt ) )+aomega 
+  +x(2,nxt)*dcos(x(6, nxt) ) /dcos( x( 4, nxt ) ) 

f (6, nxt)*x( 1, nxt ) *dsin( xC 8, nxt) )*dtan(x(4, nxt) )+x(3,nxt ) 

+  +x(2, nxt )*dcos(x(6, nxt) ) *dtan( x( 4, nxt ) ) 

return 
end 


a 


a 


a 
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PROGRAM  LISTING  torder 


subroutine  order C inoise, tseed  ) 

this  routine  reeds  the  dete  file  ’sitings’  which  was  created  by 

'truth*  and  organizes  it  in  time  order. 

the  routine  uses  a  sodified  bubble  sorting  technique 

if  noise  is  to  be  added  to  the  system,  it  is  done  here 

double  precision  tsl it ( 100 ) , sdatC 2, 100 ) 
double  precision  tper, tain, tmin2, zero, tseed 
dimension  isl it ( 100 ) , idl ( 100 ) , id2( 100 ) , iseq ( 100 ) 
constants 

tper*12. 00d+00*2. 00d+00 
open  files  and  read  data 
open( 1, files ’sitings ’ ) 
open<2, f ile= ’stars ’ ) 
if ( inoise. ne. 0)  write( * , 102 ) 

102  format(2x, //, 15x, ’NOISE  ADDED  TO  TRUE  SLIT  CROSSING  TIME’) 
k=0 

1  continue 
k=k+l 

read(l, 101)  isl it< k ) , idl( k ) , id2( k ) , tsl it( k ) , 

+  (sdat(ij.k), ij*l,2) 

101  format(6x, il,9x, 13, 4x, 14, 15x,e20. 13, /, 

+  10x,e20. 13, 10x,e20. 13,  / ) 

if  inoise  not  equal  0  noise  will  be  added 
if ( inoise. ne. 0 )  then 
call  nois( tseed, tslit(k) ) 

7  continue 
endlf 

ifddl(k).  ne.999)  goto  1 
close  (1) 

end  of  data  reading 

now  data  will  be  sorted  by  time 

tmin2  =  tsl  ltd )+tper 

do  3  i*l,k-l 

Jmin3i 

tmin-tsl it( 1 ) 
iain*l 

do  2  j*l,k-l 

if(tslit(J). lt.tmin. and. tslit(J). It.tmin2)  then 
tminstsl it( j ) 
imin*J 
end  if 

2  continue 

iseq( imin)3Jmin 

tslit(imin)*tslit(imin)+2. 0d-*-00*tper 

3  continue 

now  convert  data  back  to  original  data  and  write  to  'stars’ 
do  4  j  *1 , k-1 

tsl i t C  j ) * tsl i t ( J )-2. 00d+00* tper 

4  continue 


do  5  i*l , k-1 

if ( iaaq( i ) . eq. J  )  than 

writa(2, 101)  ial i t ( i ) , idl ( i ) ,  id2( i ) , tal i t ( 1 ) , 
♦  ( adat ( ij,i),ij*l,2) 

andif 

6  continua 
6  continua 

zaro«0. 00d+00 
izero=0 
if lag-999 

write(2, 101 )  izaro, if lag, izaro, zero, zero, zero 

return 

and 


include ’ noiae ’ 


PROGRAM  LISTING  tresid 


subroutine  resid(t,xO) 

This  routine  reeds  in  initial  data 

reads  the  stars  file  (actual  sitings) 

and  gets  an  estimated  slit  crossing  time  from  tstate 

double  precision  x0(6),t0, tf lag 

double  precision  xst(6), tst, ttrue, sdat! 2 ) 

read  star  crossing  info 

open  ( 1, file* ’stars ’ ) 

10  continue 

read(l, 100) is 1 it, idl, id2, ttrue, sdat( 1 ) , sdat(2) 

100  format (6x, il , 9x, i3, 4x, 14, 15x, e20. 13, /, 

+  lOx, e20. 13, lOx, e20. 13, / ) 

if ( idl. eq. 999)  goto  999 

cal 1  tstate ( to, ttrue, xO, sdat, islit, xst, tst ) 

write  out  true  and  estimated  conditions  to  statfile 

open! 2, file* 'statfile* ) 

write(2, 101 )  islit, idl, id2, ttrue 

101  format(2x, 'Slit  *,il,’  of  star  ' , 14, ’ ( id- ' , 14, ' )  crossed  at 
+  'true  time  :',e20.13) 

write(2,102)  sdat ( 1 } , sdat (2 ) 

102  format(5x, 'star  data  ( RA, Dec) t  ’ , 2( 2x, e20. 13 ) ) 
write(2,104)  tst 

104  format! 2x, 'ESTIMATION  TIME  ’ , / , 5x, e20. 13, / ) 

if (idl. ne. 999)  goto  10 
999  continue 

write  end  of  data  flag 
if lag*0 

tf lag*0. 00d+00 

write! 2, 101 )  if lag, if lag, if lag, tf lag 

cloae(2) 

close! 1 ) 

return 


o  a 


PROGRAM  LISTING  truth 


subroutine  truth(  to,  tend,  h,  xO,  t;f  in,  xout ) 

This  routins  is  called  by  program  main 

truth  inputs  ara  initial  tiaa  (tO),  integration  tiaestep  (h) 
and  initial  stats  (z) 

truth  outputs  ara  final  tiaa  (t)  and  state  (xout) 

This  version  of  truth  records  star  sightings  for  a  12  second 
period  froa  tO.  Each  star  is  analyzed  over  this  period,  and 
outputs  are  sent  to  'sitings',  'sitings'  stores  data  as 
idl , id2, time, star  r. a., star  dec., and  true  state  at  t 
although  true  state  is  stored,  the  estiaation  algoritha 
will  not  access  this  state 

true  state  will  be  used  only  for  result  analysis 

coaaon/haa/haat, x( 12, 4) , f ( 12, 4) , arrest ( 12) , n, hash 

double  precision  tO, z(6) , t , xout(6 ) , heat, haah 

double  precision  x,f, arrest, h 

double  precision  sdat( 2 ) , sold( 2 ) , obs( 2, 3) 

double  precision  slitx.dslit, thetai , f nterp, fov 

double  precision  tper , x0(6 ) , xf lag(3 ) , xhold(2, 3 ) , tf in, tend 

diaension  ihold(2,3) 

constants 

thetai *2. 61799388d-01 
fov*2. 617993878d-02 
tper*tend 
nf legal 

open  star  data  table  for  reading 

open( 1 , f i le* ' stardat ' ) 

open(2, f ile* 'sitings ' ) 

read  stars 

do  9  1*1,6 

z( 1 ) *x0( i  ) 

9  continue 
11  continue 

read  (1,101)  idl , ( sdat( i ) , i *1 , 2 ) , id2 
101  foraat( lx, i3,2(2x,e20. 13), 9x, 14) 
write<#,#)  idl,id2 
end  of  star  reading 
setup  data  for  haming  calls 
n*6 

t*to 

haat*t-tO 
haah*h 
do  1  i*l,n 
x(i,l)*xO(i) 

1  continue 

Initialize  haming 
write( * , * ) id2, to, t , heat 
nxt*0 

call  haming(nxt) 
if(nxt.ne.O)  goto  10 
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write( *,  2) 

2  foraat(2x, 'haaing  did  not  initialize’) 
atop 

10  continue 

if ( idl . eq. 999 )  goto  999 

aake  firat  call  to  obaer 

t»haat 

ifov*0 

iaiaa*0 

call  obaer (z, t, n£ lag, adat, oba) 

20  continue 

atore  row  2  of  oba  in  aold 
do  3  J  3 1 , 2 
aold( J )*oba( J , 2) 

3  continue 

integrate  atate  equations 
call  haaing(nxt) 
write(*,«)  idl, id2, to, t , haat 
now  call  obaer  at  new  tiae 
do  4  i*l,6 
z( i )*x( i,nxt ) 

4  continue 
t*haat 

call  obaeriz, t,nf lag, adat, oba) 
check  for  a lit  croaainga 
do  7  ialit*l,2 

alitx*aold( ialit )*oba( ialit, 2) 
if(alitx. le.0.00d+00)  then 
now  check  to  aee  if  croaaing  ia  forward  looking 
if <oba( ialit, 1 ) . gt . 0 . 00d+00 )  then 

aee  if  this  ia  the  firat  alit  croaaed 
and  advance  croaaing  counter  iaiaa 
iaiaa*iaiaa+l 
if ( ifov. eq. 0 )  then 

aee  if  the  slit  ia  in  the  sensor  fov 
dal itsdcos( fov ) 

if ( ialit. eq. 2)  dal it*dcoa( f ov/dcos( thetai ) ) 
if (oba( ial it, D.ge.dalit)  ifov-1 
write<*,«)  ialit, idl, id2,t, oba( ialit, 1), da lit, slitx 
end  if 

now  we  have  a  confiraed  siting  if  ifov.gt.O 

if ( ifov. gt. 0)  then 
hold  this  siting 
ihold( ialit, 1)* ialit 
ihold( ialit, 2  )»idl 
ihold( ialit, 3)»id2 
xhold( ial it, 2 ) «adat( 1 ) 
xhold( ialit , 3 ) 3 adat (2 ) 

need  to  interpolate  other  hold  paras  baaed  on  oba(x,2) 
k«nxt-l 


> 


► 


i 


► 


i 


£ntarp=-sold( ial it ) / ( obs( ial it , 2 )-aold( ial it ) ) 
xhold< ial it, 1 )*<  t-h )+h*fnterp 
i£ov*ifov+l 
•ndif 
•ndif 
•ndif 

7  continue 
c 

c  now  if  i£ov*3  want  to  write  atar  to  aitinga 

if ( ifov. aq. 3)  then 
do  S  iw-1,2 

write(2, 102 )( ihold( iw, iww) , iww*l, 3), ( xhold( iw, iww) , iww=l, 3) 
102  foraat( lx, ’Slit  \il,’  of  atar  ’ ,  i3,  *  (  id-  ’ ,  i4, 

+  * )  croaaed  at  t= ' , a20. 13, / , 

+  5x, *RA  =  • ,e20. 13, 5x, ‘DEC«  *,e20.13,/) 

8  continue 
•ndif 

c 

c  if  we  havent  confirmed  thia  atar  continue  integrating  to  tper 

if (t. le. tper. and. ifov. ne. 3. and. imiaa. It. 2)  goto  20 
c  if  we  havent  reached  the  end  of  the  atarfile,  get  another  atar 

if (idl.ne.999)  goto  11 
c  TEMP  END  OP  PROGRAM 
tf in*hamt+tO 
do  111  iend»l,6 
xout( iend)=x( lend, nxt ) 

111  continue 
999  continue 

c  now  write  data  flag  at  end  of  file 
do  998  1*1,3 
xf 1 ag( 1 )*0.00d+00 
998  continue 
izero*0 
if lag=999 

write(2,102)  izero, iflag, izero, ( xf lag( 1 ), 1*1,3) 

cloae< 1 ) 

cloae(2) 

c  now  compute  final  true  atate  to  pa aa  back  to  main  program 

997  continue 

call  haming(nxt) 

if (hamt. It. tend)  goto  997 

t£ in*hemt+tO 

do  996  i*l, 6 

xout( i )*x( i, nxt ) 

996  continue 
return 
end 


include ’ haming ’ 
include’ obaer v2 ’ 


PROGRAM  LISTING  tstate 

subroutine  fcstete( tepoch, t, xO, sdat, islit , xst, tst ) 
c  This  routine  uses  a  Newton-Rhaphson  method  to  estimate  a  slit 

c  crossing  time 

common/dyncom/z ( 8 ) , a( 6 ) 

common /ham/ them, xham( 12, 4 ) , f ( 12, 4 ) , errestt 12) , nham, hham 
double  precision  xst(6) , x0(6) , t,sdat(2), tst, told,tnew, tol 
double  precision  tdif , h(6),a,z,obs(2,3),dgdt 
double  precision  xham, f , errest, hham, tham 
c  constants 

nham«6 
nxt*l 

tola1.00d-08 
maxit-5 
do  1  i»l,6 
z( i )=x0( i ) 

1  continue 
iter*0 
told*t 
tO=tepoch 

2  continue 
iter=iter+l 

c  get  state  from  dyno  (1  or  2,  depending  on  order) 

call  approx2( told, to ) 
c  get  h  from  obser 

call  obser (a. told, 2, sdat , obs ) 
c  recall  left  side  of  h30 

h(  1 )*0. 00d+00 
h<2)*0. 00d+00 
h ( 3 ) =0 • 00d+00 
h(4)«obs( islit, 1 ) 
h(5)sobs( islit, 2) 
h(6)-obs( islit, 3) 
c  get  xdot  from  rhs 

do  3  i=l,6 
xham( i , nxt ) =a( i ) 

3  continue 
thamstold 
call  rhs(nxt) 

c  get  gdot  from  obser 

call  obser( a, told, 3, sdat, obs) 
dgdt*obs( isl it, 2 ) 
do  4  i»l,8 

dgdt*dgdt+h( i ) *f ( i , nxt ) 

4  continue 

c  get  g  from  obser 

call  obser(a, told, 1 , sdat , obs) 
g*obs( islit, 2) 
if (dgdt . eq. 0. 00d+00 )  then 
write( * , 101 ) 

101  formate 2x, ’PROGRAM  FAILED  in  tstate  when  dgdt  went  to  0 ’ ) 


1 
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and  if 

tnew*told-g/dgdt 
tdif *daba( tnew-told ) 

c  write(*,110)  iter , tol , tnew, told, tdif 

c  110  for»*t(2x, * debug  * , / , 5x, 16, e20. 13, / , 3(2x, #20. 13) ) 
told* tnew 

if ( iter. gt. aaxit. end. tol . ge. 1. 00d-04)  then 
write( * , 100 ) 

100  for«et(2x, 'PROGRAM  FAILED  IN  TSTATE  ( tet  couldnt  converge)’) 
stop 
endif 

if (iter.gt.aaxit)  then 
iter*0 

tol*tol*l. 00d+02 
endif 

if ( tdif . gt. tol )  goto  2 

c  if  we  have  gotten  thie  far  we  have  converged  to  an  estimate 

tst*tnew 
do  5  i*l,  6 
x«t( i )*a< i ) 

5  continue 
return 
end 


inc 1 ude ’ dy no2 ’ 


PROGRAM  LISTING  starcnv.f 


program  starcnv 
fchim  program; 

raads  star  data  from  starin 
converts  data  into  usabla  format 
orders  deta  in  scan  direction 
creates  ordered  data  file  stardat 
double  precision  atara(301 ) ,  stard(301 ) ,  rah,  ram,  dras,  pi 
double  precision  deed, decs, deem, dmin 
dimension  iseq< 3, 301 ) , bright ( 301 ) 
open( 1 , file* ’starin ’ ) 
open(2, file* ’ stardat ’ ) 
open(3,file* ’starold’ ) 
open(4, file* ’starmid ’ ) 
pi *3. 141692654d+00 
data  reading  loop 
do  1  i *1,301 

read< 1, 100)  id, irah, iram, ras, idecd, ideem, idecs, bmag 
100  format ( 12x,  14,  lx,  i2,  lx,  12,  f5. 1,  lx,  13,  lx,  12,  lx,  12,  lx,  f5.  2) 
convert  right  ascension  components  to  radians 
rah* irah 

ram* iram 
dras* ras 

stara( i )*<rah+raa/6. 00d+01+dras/ ( 6. 00d+01*6. 00d+01 ) ) 

♦  *2. 0Od+OO*pi/2. 40d+01 

convert  declination  components  to  radians 

deed* idecd 

decd*dabs( deed ) 

deem* ideem 

decs* idecs 

stard< i )*(decd+decm/6. 00d+01+deca/ ( 6. 00d+01*6. 00d+01 ) ) 

+  *pi/1.80d+02 

if ( idecd. 1 t . 0 )  stard( i ) *-stard( i ) 

here  two  stars  need  special  attention  (id’s  1766  1862) 

if ( id. eq. 1766)  stard( i )*-stard( i ) 

if ( id. eq. 1862 )  stard( i ) *-stard( i ) 

iseq(l, i )*1 

iseq( 2, i  )*id 

iseq(3, i )*i 

bright ( i )*bmag 

write(3, 200)  iseq( 1,1), stare ( i ) , stard( i ) , bright( i ) , iseq( 2, i ) , 

+  iseq( 3, i ) 

1  continue 
close( 1 ) 

now  need  to  order  the  data 
stard  goes  from  -pi  to  pi 

to  simplify  ordering  process  will  cause  stard  to  go  from  0  to  4 
do  10  1*1,301 

if (stare ( 1 ) . le. pi . and. stard ( i ) . gt . 0. 00d+00 ) 

♦  stard( i )*4. 00-stard( i )/pi 

this  takes  data  In  quadrant  4  and  assigns  it  values  between  3  &  4 


if (stara( i ) . la. pi . and. stard( 1 ) . le. 0. 00d+00 ) 

+  atard( i )*-stard( i )/pi 

c  this  takaa  data  in  quadrant  1  and  aaaigna  it  to  batween  0  and  1 
if (atara( i ) . gt. pi )  atard( i )*stard( i )/pi+2. 00d+00 
c  thia  aaaigna  data  in  quadranta  2  and  3  to  valuea  batwaan  1  and  3 
c  writa(4,200)  iaaq( 1 , i ), atara( i ), atard( i ), bright ( i ), iaaq( 2, i  ) , 

c  +  iaaq(3,i) 

10  continua 

c  now  data  haa  atard  aaaignad  in  increaaing  ordar 

c  will  aort  data  according  to  atard 

do  26  1 *1.301 
Jain*l 

d*in*stard( 1 ) 

16  iain*l 

do  20  j *1,301 

if (atard( J ) . It. dain. and. atard( J ) . It. 5. 00d+00)  then 
dain*stard( J ) 
lain* j 
andif 

20  continua 

iaaq( 1, lain)* Jain 
atard (  iain  )*atard(  lain ) +6. 00d+00 
26  continua 

c  now  all  data  haa  baan  givan  a  aaquanca  nuabar 

c  naad  to  convart  atard  back  to  valid  radian  valua  and 

c  print  aaquantially  ordarad  data  into  atardat 
do  30  i*l , 301 

atard ( i )*atard< i )-6. 00d+00 

if (atara(i). la. pi. and. atard( i ) . It. 1.00d+00)  atard( i ) =-stard< i ) *pi 
if (atara< i ) . la. pi. and. atard( i ) .  ga. 3. 00d+00 ) 

+  atard ( i )*(4. 00d+00-stard( i ) )*pi 
if (atara( i ) . gt. pi )  atard< i )*(atard( i )-2. 00d+00 ) *pi 
30  continua 

do  40  i* 1,301 
do  35  k* 1,301 
if ( iaaq( 1 , k) . aq. i  ) 

+writa( 2, 200 )  iaaq( 1 , k ) , atara( k ) , atard ( k ) , br ight ( k ) , iseq ( 2, k ) , 

♦  iaaq( 3, k ) 

200  foraat( lx, 13, 2( 2x, a20. 13) , 2x, f5. 2, 2x, 14, lx, 13  ) 

36  continua 
40  continua 

c  now  wrlta  flag  at  and  of  atardat  and  cloaa  file 

iaeq< 1, 1 )*999 
atara( 1 )*0. 00d+00 
atard ( 1 ) *0 . 00d+00 
bright ( 1 ) *0. 00 
laaq(2, 1 )*0 
laeq(3, 1 )*0 

writa( 2, 200 ) iaaq( 1,1),  atara( 1 ) , atard ( 1 ) , bright( 1 ) , iseq (2, 1 ) , 

♦  iaeq( 3, 1 ) 
cloaa(2) 
and 


TEMPORARY  FILE 


sitings 


Slit  1  of  star  77  (  id-2261)  croasad  at  t  =  0. 9616549653323e+00 
RA  *  0. 1617505701057e+01  DEC=  -0. 1304604527215e+01 

Slit  2  of  star  77 ( id-2261)  crossed  at  t*  0. 9742219666100e+00 
RA  *  0. 1617605701067e+01  DEC*  -0. 1304604527215a+01 

Slit  1  of  star  130( id-6746)  crossed  at  t*  0. 3453228588421e+01 
RA  *  0. 4733383837519e+01  DEC*  -0. 5310212731254e+00 

Slit  2  of  star  130( id-6746)  crossed  at  t*  0. 3468579114165e+01 
RA  *  0. 4733383837519e+01  DEC*  -0. 5310212731254e-K>0 

Slit  1  of  star  133( id-6742)  crossed  at  t*  0. 3481622288535e+01 
RA  *  0. 4729973173213e+01  DEC*  -0. 5162974816282e+00 

Slit  2  of  star  133( id-6742)  crossed  at  t*  0. 3496661445198e+01 
RA  *  0. 4729973173213e+01  DEC*  -0. 5162974816282e+00 

Slit  1  of  star  160( id-6752)  crossed  at  t*  0. 4552968869622e+01 
RA  *  0. 4732772972169e+01  DEC*  0. 4367201640005e-01 

Slit  2  of  star  160( id-6752)  crossed  at  t*  0. 457294460948Oe+Ol 
RA  *  0. 4732772972 189e-*-01  DEC*  0. 4367201640005e-01 

Slit  1  of  star  168( id-6771)  crossed  at  t*  0. 4788964828677e+01 
RA  *  0. 4741252363620e+01  DEC*  0.  1668680209229e-*-00 

Slit  2  of  star  168( id-6771)  crossed  at  t*  0. 4806517806391e+01 
RA  *  0. 4741252363620e+01  DEC*  0. 1668680209229e+00 

Slit  1  of  star  177( id-6787)  crossed  at  t=  0. 51651 11553339e+01 
RA  =  0.4747710081728e+01  DEC*  0. 3632272580715e+00 

Slit  2  of  star  177( id-6787)  crossed  at  t=  0. 5179998187776e+01 
RA  =  0. 4747710081728e+01  DEC*  0 .  3632272580715e-*-00 

Slit  1  of  star  219C id-6927)  crossed  at  t*  0. 6898057723155e+01 
RA  =  0. 4805487752251e+01  DEC*  0.  1269314939363e-*-01 

Slit  2  of  star  219( id-6927)  crossed  at  t=  0. 6911541012503e+01 
RA  *  0. 4805487752251e+01  DEC*  0. 1269314939363e+01 

Slit  1  of  star  232<id-  424)  crossed  at  t=  0. 7484197884076e+01 
RA  *  0.  5930192465628e+00  DEC*  0. 1556731882109e+01 

Slit  2  of  star  232(id-  424)  crossed  at  t=  0. 7502045042247e-*-01 
RA  *  0. 5930192465628e+00  DEC*  0. 1556731882109e+01 

Slit  1  of  star  243( id-2165)  crossed  at  t=  0. 8273669495085e+01 
RA  *  0. 1620065517349e+01  DEC*  0.  1 147083714065e+01 


> 

I 


r 


Slit  2  of  star  243( id-2165)  crossed  at  t*  0. 83016S8652403e-*-01 
RA  *  0. I620065517349e+01  DEC*  0. 1147083714065e+01 

Slit  1  of  star  248( id-2077)  crossed  at  t=  0. 8658927419005e+01 
RA  *  0. 1563160511522e+01  DEC*  0. 9474471364 320e+00 

Slit  2  of  star  248( id-2077)  crossed  at  t*  0. 8671526066400e+01 
RA  *  0. 1563160511522e+01  DEC*  0. 9474471364320e+00 


1 


' 


► 

D  2 


/  *,»  V  V  V 


TEMPORARY  FILE 


stars 


1  77  2261 

0. 16 1750570 1057e+01 

2  77  2261 

0. 1617505701067e+01 

1  86  7228 

O. 54726107700516+01 

2  86  7228 

0. 54726107700516+01 

1  130  6746 

0.47333838375196+01 

2  130  6746 

O. 4733383837519e+01 

1  133  6742 

0.47299731732136+01 

2  133  6742 

O. 47299731732136+01 

1  160  6752 

O. 4732772972169e+01 

2  160  6752 

O. 4732772972169e+01 


0. 9623484884 174e+00 
-0. 1304604527215e+01 

0. 9767395951997e+00 
-0. 13046045272 15e+01 

O. 14804967858736+01 
-0. 1553658 163371 e+01 

0. 1 49575960367 le+01 
-O. 1553658 16337 le+01 

0. 3453947896009e+01 
-0. 53102127312546+00 

0. 347033685004 le+01 
-0. 5310212731254e+00 

0. 3482325372939e+01 
-O. 5162874816282e+00 

0. 3500374983903e+01 
-O. 5 1629748 16282e+00 

0. 4553633795 184e+01 
O. 4367201640005e-01 

0. 457381 1172263e+01 
0. 4367201640005e-01 


168  6771 

O. 4741252363620e+01 


168  6771 

0. 4741252363620e+01 


177  6787 

0. 4747710081728e+01 


177  6787 

O. 4747710081728e+01 


219  6927 

0. 480548775225 le+01 


219  6927 

0. 48054877522516+01 


232  424 

0. 6930192465628e+00 


O. 4789667475850e+01 
0. 1668680209229e+00 


0. 4806009442064e+01 
0. 1668680209229e+00 


0. 5165846987642e+01 
0. 3632272580715e+00 


O. 5180163123040e+01 
0. 3632272580715e+00 


0.68987859436186+01 
0. 1269314939363e+01 


0. 6910557671146e+01 
0. 1269314939363e+01 


0.74849067158386+01 
0. 1 55673 1 882 109e+01 


D3 


232  424 

0. 59301924656286+00 

243  2165 

0. 1620065617349e+01 

243  2165 

0. 1620065517349e+01 

248  2077 

0. 15631605115226+01 

248  2077 

0. 1563160511522e+01 

255  2088 

0. 15637786489650+01 

255  2088 

0. 15637786489650+01 

270  2134 

0. 15846644223360+01 

270  2134 

0. 15S4664422336a+01 

283  2169 

0. 15999724143880+01 

283  2159 

O. 15999724143880+01 

999  0 

0.  o+OO 


0. 76009335554330+01 
0. 15567318821090+01 

0.82744307602160+01 
0. 11470837140650+01 

0. 83007512044910+01 
0. 11470837140650+01 

0.86696036691530+01 
0. 9474471364320O+00 

0.86706674724660+01 
0. 9474471364320e+00 

0. 89710523135390+01 
0.78447216936900+00 

0. 8979626297276O+01 
0.78447216936900+00 

0. 96933204503130+01 
0.40605569866630+00 

0. 97063275924320+01 
0. 40605669866630+00 

0. 99764228639100+01 
0. 25780462310050+00 

0. 99959050901 70e+01 
0.25780452310050+00 

0.  o+OO 

0.  o+OO 


TEMPORARY  FILE  statfile 

Slit  1  of  star  77( id-2261)  crossed  at  true  time  :  0. 9623484884174e+0 
star  data  ( RA,  Dec  )  t  0. 1617506701057e+01  -0.  1304604527216e«-01 

ESTIMATION  TIME 

0. 9623493674646e+00 

Slit  2  of  star  77( id-2261)  crossed  at  true  time  :  0. 9767395961997e+0 
star  data  (RA, Dec) t  0. 1617505701057e+01  -0. 1304604627216e*01 

ESTIMATION  TIME 

0. 9767361246182e+00 

Slit  1  of  star  86( id-7228)  crossed  at  true  time  t  0. 1490496788873e+0 
star  data  (RA.Dec)t  0. 5472610770051e-*-01  -0.  1553658163371e+01 

ESTIMATION  TIME 

0. 149049315800Se+01 

Slit  2  of  star  86(  id-7228)  crossed  at  true  time  :  0. 149675960367le+0 
star  data  (RA,Dec)>  0. 6472610770051e+01  -0.  1563668163371e-«-01 

ESTIMATION  TIME 

0. 1495756293484e+01 

Slit  1  of  star  130( id-6746)  crossed  at  true  time  *  0. 3453947896009e+0 
star  data  ( RA, Dec ) i  0. 4733383837519e+01  -0 . 5310212731254e+00 

ESTIMATION  TIME 

0. 346395331 4 138e+01 

Slit  2  of  star  130( id-6746)  crossed  at  true  time  i  0. 347033685004 le+0 
star  data  (RA,Dec)«  0. 4733383837619e+01  -0. 5310212731264e+00 

ESTIMATION  TIME 

0. 3470366689066e+01 

Slit  1  of  star  133( id-6742 )  crossed  at  true  time  :  0. 3482325372939e+0 
star  data  ( RA, Dec ) s  0. 4729973173213e+01  -0. 5162974816282e+00 

ESTIMATION  TIME 

0. 348233241 8755e+01 

Slit  2  of  star  133( id-6742)  crossed  at  true  time  :  0. 3500374983903e+0 
star  data  ( RA , Dec ) :  0. 4729973173213e+01  -0. 5162974816282e+00 

ESTIMATION  TIME 

0. 3500363708282e+01 

Slit  1  of  star  160<  id-6762)  crossed  at  true  time  :  0. 4563633795184e-*-0 
star  data  ( RA, Dec ) :  0. 4732772972169e+01  0. 4367201640005e-01 

ESTIMATION  TIME 

0. 456362698e230e+01 

Slit  2  of  star  160( id-6752)  crossed  at  true  time  :  0. 467361 1172263e+0 
star  data  (RA,Dec)i  0. 4732772972 169e+01  0. 4367201640005e-01 

ESTIMATION  TIME 

0. 4573609821861e+01 


Slit  1  of  star  168( id-6771)  crossed  at  true  time  :  0. 4789667475850e+0 
star  data  (RA,Dec):  0. 474l252383620e+01  0. 1668680209229e+00 

ESTIMATION  TIME 

0. 4789669653552e+01 

Slit  2  of  star  168( id-6771)  crossed  at  true  time  :  0. 4806009442064e+0 
star  data  (RA, Dec):  0. 474l252363620e+01  0. 1668680209228e+00 

ESTIMATION  TIME 

0. 4805987380 158e+01 

Slit  1  of  star  177( id-6787)  crossed  at  true  time  s  0. 6166846987642e+0 
star  data  ( RA, Dec ) :  0. 4747710081728e+01  0. 36322725607 16e+00 

ESTIMATION  TIME 

0.  5185861588250e+01 

Slit  2  of  star  177(  id-6787)  crossed  at  true  time  :  0. 5180l63123040e-*-0 
star  data  ( RA, Dec ) :  0. 4747710081728e+01  0. 3632272680715a+00 

ESTIMATION  TIME 

0. 5180175489028e+01 

Slit  1  of  star  219( id-6927)  crossed  at  true  time  *  0. 6898785943618e+0 
star  data  ( RA, Dec ) *  0. 4806487762251e+01  0.12693149393630+01 

ESTIMATION  TIME 

0. 6898787988382e+01 

Slit  2  of  star  219( id-6927)  crossed  at  true  time  t  0. 6910567671 146e+0 
star  data  (RA,Dec)»  0. 480548776225 le+01  0. 12693149393630+01 

ESTIMATION  TIME 

0. 69 1066021 0426e+01 

Slit  1  of  star  232(id-  424)  crossed  at  true  time  :  0. 74849067 15838e+0 
star  data  ( RA, Dec ) :  0. 5930192466628e+00  0. 1566731882109e+01 

ESTIMATION  TIME 

0. 7484890 136374e+01 

Slit  2  of  star  232(id-  424)  crossed  at  true  time  :  0. 7500933555433e+0 
star  data  ( RA»  Dec ) :  0. 5930192465628e+00  0. 1566731882109e+01 

ESTIMATION  TIME 

0.  75009 12603769e-*-01 

Slit  1  of  star  243( id-2165)  crossed  at  true  time  :  0. 8274430760216e+0 
star  data  ( RA , Dec ) :  0. 1620065517349e+01  0. 1147083714065e+01 

ESTIMATION  TIME 

0. 8274426586969e+01 

Slit  2  of  star  243( id-2165)  crossed  at  true  time  :  0. 8300751204491e+0 
star  data  (RA,Dec)>  0. 1620066517349e+01  0 . 1 147083714065e+01 

ESTIMATION  TIME 

0. 83007485702 14e+01 


Slit  1  of  star  248( id-2077)  crossed  at  true  time  :  0. 8659603669153e+0 
star  data  ( RA, Dec ) t  0. 1663160611622e+01  0. 9474471364320e+00 


ESTIMATION  TIME 

0. 86596 15809687a+01 

Slit  2  of  star  248 (id-2077)  crossed  at  true  time  *  0. 8670657472465e+0 
star  data  (RA, Dec) t  0. 166316061 1522e+01  0 . 947447 1364320*+ 00 

ESTIMATION  TIME 

0. 8670666897236*+01 

Slit  1  of  star  255( id-2088)  crossed  at  true  time  :  0. 897l052313539e+0 
star  data  (RA,Dec)*  0. 1563778648966e+01  0. 7844721693690e+00 

ESTIMATION  TIME 

0. 897107 1742029e+01 

Slit  2  of  star  265( id-2088)  crossed  at  true  time  :  0. 8979626297276e+0 
star  data  (RA, Dec):  0. 1563778648965e+01  0. 7844721693690e+00 

ESTIMATION  TIME 

0. 8979646825753e+01 

Slit  1  of  star  270( id-2134)  crossed  at  true  time  :  0. 96933204603 13e+0 
star  data  ( RA, Dec ) :  0. 1684684422336e+01  0 . 4060666986663e+00 

ESTIMATION  TIME 

0. 9693318591079e+01 

Slit  2  of  star  270( id-2134)  crossed  at  true  time  :  0. 9706327592432e+0 
star  data  ( RA , Dec ) :  0. 1584864422336e+01  0 . 4060653986663e+00 

ESTIMATION  TIME 

0. 9706322562576e+01 

Slit  1  of  star  283( id-2169)  crossed  at  true  time  »  0. 99764228639 10e+0 
star  data  ( RA, Dec ) :  0. 1599972414388e+01  0. 2578045231005e+00 

ESTIMATION  TIME 

0. 99764 1233 1844e+01 

Slit  2  of  star  283( id-2159)  crossed  at  true  time  :  0. 9995905090170e+0 
star  data  (RA, Dec) :  0. 1599972414388e+01  0. 2578046231006e+00 

ESTIMATION  TIME 

0, 999689553076 le+01 


Slit  0  of  star 


0(  id- 


Ol  crossed  at  true  time  :  0. 


e+0 


TEMPORARY  FILE 


dtfilc 


Slit  1  0. 9623484884174«+00  0. 962340387464e«+OO 

0.4884442720068«-01  -0. 10322O6169043«-O1  0. 1842105428988®+01 

0. 81000 18988367«-01  -0. 6827933647731«-01  0. 1911392043842«+01 

Slit  2  0. 978730696 1997«+00  0. 9767361246182«+00 

0. 45698 18340783«-01  -0. 5121896365819«+00  0. 1903776241063«+01 

0. 5910863206373*+00  -0. 1973604216770«+OO  0. 1840672941823«+01 

Slit  1  0. 1490496786873«+01  0. 1400493 168008«+01 

0. 1657492248603«+00  -0. 2786672914790«-01  0. 28473 11 453673«+01 

0.9201886675446«-01  0. 2424186693488«-01  0. 1907344166192«+01 

Slit  2  0. 1496769803871«+01  0. 1496755293484«+01 

0 .  1573668682420«-*-00  -0 . 7392385874678«+00  0 .  2908219073380«+01 

0. 6208710671477«+00  0. 3098121074822«-01  0. 10294884396O7e+Ol 

Slit  1  0.3463947898009«+01  0. 3453953314138«+01 

0. 3798249834843«+00  0. 134931 1647171e+00  0. 86831 08660669«+01 

-0. 466075238367 1«-01  0. 6487730619669«-01  0. 190863 1397698a+01 

Slit  2  0. 3470336850041e+01  0. 3470368589056«+01 

0 .  3777338689896«+00  -0 .  8280423090686«+00  0 .  6680762625660«f  0 1 

0. 2169407732724«+00  0. 5126096866867«+00  0. 18938 14760742«+01 

Slit  1  0.3482325372939«+01  0. 3482332418765«+01 

0.  371 287254 1496«+00  0. 1393 1011 12686«+00  0. 6837783825328«+01 

-0. 428062022480 la-01  0. 8083856757213«-01  0. 1906680830 178«+01 

Slit  2  0.3600374983903«+01  0. 3500383708282«+01 

0 . 3882573579883«+00  -0 . 8 167263368482«+00  0 . 6738283320765«+0 1 

0. 2118472874077«+00  0.  51 163 14532844**00  0.  1893246670447*-*-01 

Slit  1  0.4553633795184*+01  0. 4683628988230«+01 

0 . 4098308389756*+00  0 . 3014775059426**00  0 . 8676865097808«+01 

-0. 8190872135005*— 01  -0. 23080484 15446*-01  0. 1911244868251**01 

Slit  2  0.4573611172263*+01  0. 4873809821861«+01 

0.4049838617559*+00  -0 . 3620724 168606**00  0. 8761189830986*+01 

-0.  104141 1707509*+00  0.  4857225785489**00  O. 1882869397685*+01 

Slit  1  0.4789667475860*+01  0 . 4789869853652«^01 

0. 42538822948 14*+00  0. 3296067896979*+00  O. 9126175587601*+O1 

-0.9142098152716*-01  -0. 4733420035519*-01  0. 1913162384894**01 

Slit  2  0. 4806009442064**01  0. 4805987380158**01 

0.4211287762175*+00  -0. 2434168931074**00  0. 9197154369800*+01 

-O.  1754 135069 165*+00  0 . 4640645364879** 00  0. 1882601787411*+01 

Slit  1  0.  5165846967642O01  0. 6165861888250e+01 

0.4222631 33839 l«+00  0. 36406021 13571**00  0. 9844020187381**01 

-0.8080188807740*-01  -0. 8213376281860*-01  0.1915120156626**01 

Slit  2  0. 6180163123040**01  0. 5180175489028**01 

0.4178252646043*+00  -0. 4464412348864*-01  0.98-  38661790**01 

-0.2600290197687**00  0. 39230348791 11**00  0. 1882ui6077163*+01 

Slit  1  0. 6898785943618**01  0. 6898787988382e*01 

0.4171670199318**00  0. 3773414278627*+00  0.1316518757451**02 

0. 6918974309387*-01  -0 . 9024440145628*-01  0. 1913269772024*+01 

Slit  2  0. 6910657671146**01  0. 6910560210426*+01 

0.41 335828608 16*+00  0. 8085862535571*+00  0. 1315835468950*+02 

-0.4105338359647*+00  0 . 5862147960830e-01  0. 1884909229074**01 

Slit  1  0.7484906716838*+01  0. 7484890136374**01 

0.4941750996559**00  0. 367362l424956*+00  0. 1428366558302**02 


i 


TEMPORARY  FILE 


0. 

0. 

0. 

0. 

0. 

-12.0000000 
-24. 0000000 
-36.0000000 
-48.0000000 
-60.0000000 
-72.0000000 
-84.0000000 
-96.0000000 
108.0000000 
120.0000000 
132.0000000 
144.0000000 
156.0000000 
168.0000000 
180.0000000 
192.0000000 
204.0000000 
216.0000000 
228.0000000 
240. 0000000 
252.0000000 
264. 0000000 
276. 0000000 
288.0000000 
300. 0000000 
312.0000000 
324.0000000 
336. 0000000 
348.0000000 
360.0000000 
360.0000000 
360. 0000000 
360.0000000 
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